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This paper is an English translation of Michel Henon's thesis, Sur revolution 
dynamique des amas globulaires originally published in French in the Annales 
d'Astrophysique, Vol. 24, p.369 (1961). 

Conventions and notations are as in the original version, for consistency. 
The English version is written so that it is as faithful to the French text as 
possible. The translator added some notes [T.N.] for the sake of clarity, when 
required. Original French, English and Russian abstracts written by M. Henon 
are available in the original version of the paper. The English part is reproduced 
below. 
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As the author explains himself on the cover page of the original version: 
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de Paris. The thesis has been defended in Paris, on 1961 December 11, before 
a jury composed of Messieurs Danjon (president of the jury), Schatzman and 
Delcroix. This work is not related to a series of papers previously published 
and entitled L 'amas isochrone. 



ORIGINAL ABSTRACT 

Chapter I: The structure and evolution of globular clus- 
ters are closely linked; we propose here to study them si- 
multaneously, with a purely theoretical approach. The es- 
sential hypotheses are: (1) spherical symmetry; (2) quasi- 
permanent regime; (3) isotropy of the velocities at all 
points; (4) mass equality. 

Chapter II: We establish the system of fundamen- 
tals equations (2.25). The cluster model is reduced to 
a canonic form by a homology transformation with time- 
dependent parameters. We will look for a model of invari- 
able canonic form, that is a model which remains similar 
to itself while evolving. 

Chapter III: We show that the galactic field imposes 
the relation (3.8): the cluster's radius is proportional to 
the cube root of its mass. 

Chapter IV: Preliminary calculations show that the 
density must be infinite at the center of the cluster. We 
obtain the asymptotic expression (4.14) for the distribu- 
tion function near the center. The conservation of mass 
imposes condition (4.32). The existence of an energy flux 
toward the cluster center is predicted. 

Chapter V: The ensemble of equations and conditions 
obtained in the preceding chapters forms a system which 



is resolved numerically by means of an electronic calcu- 
lator. We find a unique solution: the "homology model" 
[T.N.: homologous model] (Table ffl Figures H El [U ED 
I16p . Its mass and radius are finite. The outer radius is 
approximately 10 times the mean radius. The mass de- 
creases linearly as a function of time. About one-third 
of the negative energy of the cluster is carried off by the 
stars which escape; the other two-thirds accumulate in 
the center, in multiple stars. 

Chapter VI: We introduce a small number of stars of 
different mass in the homology model and we calculate 
their distribution (Figures [8] and |9]) and their escape rate 
(Figure [TU)) . We find that as their mass increases, they 
are more concentrated, and escape less rapidly. In par- 
ticular, if their mass is greater than 3/2 the mass of the 
normal stars, they are nearly all collected near the center 
of the cluster, and their escape rate is zero. 

Chapter VII: We show that should the initial central 
density of a cluster be finite, it will grow and become in- 
finite within a finite time (Figure [T2")) . We study next the 
evolution of a cluster differing slightly from the homology 
model: we find that the differences decreases, whatever 
their form. Thus the homology model is very probably 
the final state toward which the clusters tend. 
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Chapter VIII: The theoretical results are compared 
with the observational data on globular clusters. First 
we construct an "artificial cluster" (Figures [15] and [T6|) 
which permits a general comparison. Then we compare 
the projected densities in detail, Sandage's star counts 
for the cluster M3, and the brightness measurements for 
47 Tuc by Gascoigne & Burr. The agreement is very sat- 
isfying (Figures HUGE ED- For the stars of M3 > we nnd a 
definite mass-luminosity relation, which however is in dis- 
agreement with the relation indicated by stellar evolution 
theory (Figure 1201) . The observed masses and radii are in 
good agreement with the theoretical relation (Figure l2"2l . 
The outer radii of the clusters are about twice as large 
as the observed radii. The absolute escape rate is con- 
stant and equals: 2.3 x 10 -6 M Q / yr. This constancy is 
confirmed by the observed distribution of globular cluster 
masses fFigure l2"3"]l . We obtain the initial mass distribu- 
tion. Observation shows that uj Cen is the only globular 
cluster whose central density is not yet infinite. From this 
we deduce that the age of the globular clusters must lie 
between 24 x 10 9 yr and 44 x 10 9 yr, in good agreemen 
[sic] with the age deduced from stellar evolution. 

Chapter IX: We list a number of directions in which 
development of research is to be hoped for. 

1 INTRODUCTION 

The problem of the dynamical evolution of a globular clus- 
ter can be stated very simply. The only important force 
is the mutual attraction between the stars of the cluster; 
the other forces (like radiative pressure, electromagnetic 
forces, relativistic effects, etc.) are negligible. There- 
fore, the topic is the classical n-body problem: finding 
the motion of n points of given masses, mutually attract- 
ing themselves as the inverse square of their distance. 

This exposition, whereas simple, relates to an ex- 
tremely arduous mathematical problem. Despite a large 
number of studies, it has not been possible to find an ex- 
plicit solution, which very likely does not exist. Hence, 
one can think of t he numerical integration. This way, 
von Hoerner (1960) computed the evolution of artificial 
clusters comprising up to 16 stars, thanks to an electronic 
device. But high values of n are out of reach with such a 
method, as the computational time becomes rapidly ex- 
treme, even for a machine; the case of n = 16 already 
corresponds to a system of 96 simultaneous differential 
equations. 

In the globular clusters, n is of the order of magnitude 
of 10 5-6 . Such a high value naturally suggests to give up 
following the individual motions, and to use a statistical 
method. The structure of the cluster is then defined at all 
times thanks to a distribution function in a 7-dimensional 
space: 3 position coordinates, 3 velocity coordinates, and 
the mass. One can write a system of equations that al- 
lows, in principle, to calculate the evolution of this func- 
tion, given its initial form. Unfortunately, the numerical 
resolution seems absolutely out of reach with this gen- 
eral form: indeed, one faces an integro-differential system 



with 8 independent variables! 

Thus, it is necessary to make simplifications, or ad- 
ditional hypotheses. Two of them are classical and well 
valid: 

(HI) we assume that the cluster is spherically symmetric; 
(H2) we assume that the cluster has reached the steady 
state: 

The number of independent variables is then reduced from 
8 to 4 (see Equation 2.2 below); but the equations are still 
too involved. Therefore, one must make new simplifica- 
tions, much more arbitrary, that are often only justified 
by their practical utility. The main one consists in consid- 
ering the structure and the evolutio n problems separately . 
This is how a ser i es of studies (e.g. |Plummerl[l~91ll Il915l : 
lEddingt onl Il916t iJeansI 119 16: Chandrasekhar 1942 Sec- 
tion 5.8. ICamra|[l952t IWoollev fc RobertsoiJll956t iHenonl 
19591 ) has focussed on the structure of the clusters, with- 
out considering the evolution. In this case, the distri- 
bution function can take any form; it is set by an addi- 
tional, arbitrary hypothesis, changing from one author to 
the other. The models obtained this way represent well 
the observed clusters; however this agreement does not 
mean much, because one of the fundamental equation of 
the problem as been suppressed, and replaced by an ad 
hoc hypothesis. 

Other works have, on the contrary, studied the evo- 
lution by assumin g a known structure for the clust er 
among oth ers, see lSnitzerlll940|: Chandrasekhar 1943b. c; 



1956 



KinJl958allb1|cllSpitzer fc Harmlll958tlvon HoernerL 
Agekiarj|l959t lHenoiJll96d: lKingj|l960l: lMichieJll96l[ ). In 

most of the cases, it is assumed that the cluster is homo- 
geneous and that the gravitational potential is uniform; 
the advantage of this is to get rid completely of the space 
variables. However, the results obtained (the major one 
being the escape rate of the stars from the cluster) often 
disagree from one author to the next (see Figure [TU|) , as 
a consequence of the arbitrary hypothesis, that has been 
introduced here again. 

In fact, the two aspects, structure and evolution, can- 
not be told apart; they are intimately related. Equations 
show this clearly (see Equation 2.22 below). First solv- 
ing the problem of the structure of the cluster, and then 
finding out its evolution, or the other way round, is not 
possible. All the equations have to be treated and solved 
as a whole. 

We hope that the present work constitutes a first step 
toward this. The homologous model presented in Chap- 
ter V is obtained by solving simultaneously the equations 
of structure and evolution. However, this model is far 
from being a satisfactory solution to the problem. Indeed, 
it has been necessary to keep two arbitrary hypotheses, 
to avoid a too high complexity in the calculations: 
(H3) we assume that the distribution function only de- 
pends on the total energy (this is equivalent to assuming 
that the velocity distribution is everywhere isotropic); 
(H4) we assume that all the stars have the same mass. 

It is difficult to known how these simplifications affect 
the correctness of the results. Therefore, the conclusions 
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should be taken carefully, and considered as illustrations 
instead of definitive statements. In fact, the main goal 
of this work is not to immediately obtain useful informa- 
tions, but rather to set a more rigorous theory by getting 
rid of one of the arbitrary postulates it used to have, 
and by pointing at those that remain to be narrowed (see 
Chapter IX). 

One can ask whether it would be possible to make 
progress thanks to direct observations, and thus to cir- 
cumvent or at least diminish, the difficulties of the theo- 
retical work. Unfortunately, it seems that such a method 
is not very successful. Observing globular clusters is dif- 
ficult, because of a number of reasons: the stellar density 
varies a lot from the center to the outskirts, so that it is 
impossible to get a satisfactory image of all the parts of 
the cluster on a single plate; one only observes the pro- 
jected density and not the spatial density; measurements 
of individual proper motions are impossible to make; mea- 
surements of radial velocities are very imprecise, because 
of large distances. Worst, one observes only the brightest 
stars of the cluster; these stars only contribute to a small 
fraction of the total mass. In other words, the essential 
structure of the cluster remains invisible. 

Observations of globular clusters bring too little infor- 
mation if one is to derive the theoretical model, charac- 
terized by a multi-variable function. This point can be 
illustrated thanks to the following quotes about the com- 
parison between theoretical models and observed clus- 
ters: [T.N.; in english in the original version] "the law 

represents the structure of a globular cluster with 
close approach to accuracy everywhere, with the excep- 
tion of the regi on immediately surrounding the center" 
( Plumme"rlll911 ); "this model appear s to fit the var iations 
of brightness of the observed disk" (|Cammlll952l ); "cal- 
culations with this model give projected densities which 
agree much better with ob servations than do those o f 
the simple isothermal case" (|Woollev fc Robertsorll956h : 
"with the first mo dels we can obtain a rather fair fit" 
( von Hoernerlll957 ): "l'amas isochrone ... ressemble d'une 
maniere surprenante aux amas reels" [T.N.: the isochrone 
cluster ... loo ks surprisingly similar to the real clusters] 
( Henonl 1959 ). But each time, the model is different... 

Thus, we can say that the agreement with observa- 
tions is a necessary, but not sufficient condition for the 
validity of a theory on the dynamics of the clusters. Ob- 
servations can be a final test for the theory, but not a 
starting point. 

However, the theory can precisely progress in a purely 
deductive way, only owing to its own strengths. As we 
have seen, the physical problem can be translated into a 
system of equations describing it, if not perfectly rigor- 
ously, at least with a very good precision. This system 
encloses all the information required to solve the prob- 
lem. That is, we feel that efforts should be made in this 
direction. Obviously difficulties exist; but they are only 
technical, mathematical difficulties, not conceptual ones. 
Thanks to the powerful computing machines nowadays, it 
does not seem impossible to solve the problem anymore. 



2 EQUATIONS 

2.1 Definitions and hypotheses 

Let 

<p(x, y, z, v x , v y ,v z , m, t) dx dy dz dv x dv y dv z dm (2.1) 

be, at time t, the number of stars of mass between m and 
m+dm and whose the six spacial and velocity coordinates 
lie between x and x + dx ... v z and v z +dv z . As we will see, 
this distribution function (p fully describes the structure 
of the cluster. 

We assume that the cluster is spherically symmet- 
ric and that it has reached a steady state or, strictly 
speaking, a quasi-ste ady state because evolution makes 

"~ 19401 : iKuzminl 



it change slowly (see IS pitzcr 



can show (Jeans 1915; Kurth 



15; 



1957). One 



1955) that the distribution 



function depends on the position and velocity coordinates 
only through two invariants (or strictly speaking, quasi- 
invariant): 

A, the angular momentum of the star; 
E, the total energy of the star, per unit mass (we will 
refer to this quantity as total energy, for conciseness). 
Therefore, we have: 



tp(x, y, 2, v x , v y , v z ,m, t) = f(E, A, m, t). 
E is given by the fundamental relation: 



E 



v 



(2.2) 



(2.3) 



where U is the gravitational potential, and v is the veloc- 
ity of the star. 

Here comes our first arbitrary hypothesis. We assume, 
for the rest of the paper, that / does not depend on the 
angular momentum A. This is equivalent to assuming 
that the velocity distribution is isotropic (i.e. spherically 
symmetric) everywhere in the cluster, as shown in (2.3). 
Thus, the distribution function shrinks to: 



f(E,m,t). 



(2.4) 



Thanks to this hypothesis, one can have the following 
reasoning. The mutual perturbations of the stars tend to 
establish a Maxwellian, therefore isotropic, velocity dis- 
tribution everywhere. The escape phenomenon goes, as 
we know, against such a distribution: every star rising 
above a certain critical velocity leaves the cluster (see 
Chapter III). But this limitation only concern the modu- 
lus of the velocity and not its orientation; in other words, 
it is itself spherically symmetric. We can therefore as- 
sume that the Maxwellian equilibrium, which is impossi- 
ble to establish in term of the distribution of the velocity 
moduli, does occur in term of the distribution of the ori- 
entations, i.e. the velocity distribution becomes isotropic 
everywhere, after a sufficient time. 

But this reasoning is ambiguous: the stars travel 
across the entire cluster, and thus, it is forbidden to con- 
sider the evolution of a small region in isolation. There- 
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fore, the validity of our hypothesis is not ensured!]]- We 
will examine the perspectives about this point in Chap- 
ter IX. 

2.2 Equations of structure 

These equations are well-known. The density at a given 
point comes from the distribution function through: 



mip dv x dv y dv z dm, 



or, using the isotropy of the velocities and (2.3): 



(2.5) 



p = 4n mdm (2E - 2U) 1/2 fdE. (2.6) 
Jo Ju 

Furthermore, p, U and the distance to the center r are 
linked through Poisson's equation: 



d 2 u 

dr 2 



2dU 
r dr 



4irGp. 



(2.7) 



We note that U depends on both r and the time t because 
of the evolution of the cluster: U — U(r,t). In the same 
way: p = p(r, t). The equations (2.6) and (2.7) give p and 
U, i.e. the structure of the cluster, once the distribution 
function / is known. It remains to write an essential 
equation: the one ruling the variation of / with time. 

2.3 Local evolution equations 

If the potential of the cluster was perfectly regular, every 
star would keep its total energy forever and no evolution 
would occur. In reality however, the potential is created 
by a discrete distribution of masses and yields irregular 
features, that are random and change with time. As a 
consequence, the total energy of the star experiences per- 
tur bations. This question has often been addressed (see 
e.g. IChandrasekharlll942L Section 2.1); one showed that 
the star mostly experiences many small perturbations, 
whose cumulative effect leads to a slow, quasi-continuous, 
variation of the total energy of the star. 

This total effect is proportional to ln(//Zi), where I is 
the maximum distance to the perturber stars, and l\ is 
the "impact parameter" , i.e. the distance by which two 
stars have to be separated to be deviated by an angle 
of 90°. That is, even the very distant stars have a non- 
negligible effect. However, l\ is very small: one can show 
that: 

h * ~, (2.8) 
n 

where r e is the radius of the cluster. (Recall that n is the 
total number of stars.) Consequently, the majority of the 
perturbations are caused by relatively close stars: e.g., if 
the cluster counts 10 5 stars, l\ equals r e x 10~ 5 , and 80% 
of the perturbations would come from stars closer than 
r e x 1CP 1 . Therefore, in order to calculate the perturba- 
tions at a given point in the cluster, we will admit that 



iMichid dl96lT) highlights an increase of the anisotropy in the 
external regions, for a particular model. 



the distribution of the velocities is the same in this point 
as everywhere else. This approximation probably leads 
to an error of about 10%, which is acceptable given the 
current state of the theory. 

The distribution function depends on r and v, through 
E; to emphasize this, we set: 

v 
~2 



f{E,m,t) = f [U+—,m,t\ =a(r,v,m,t). (2.9) 



In this entire Section, we will study what happens 
at a given point in the cluster, at the distance r from 
the center. The function a defined in (2.9) describes the 
vel ocity distribution at this point. 

Rosenbluth. MacDonald fc Juddl (|l957h gave the 



equation of evolution of such a distribution, with the ef- 
fects of perturbations, in the very general case where there 
is no symmetry for the velocities, as well as in the case 
of axisymmetry. In our case of spherical symmetry for 
the velocities, their equation becomes (after some calcu- 
lations jfp 



da \ 

to), 



= 16tt 2 G 2 

]_d_ 

v 2 dv 
mi da 



In (n) / midrai 
'o 



/ a\v\dvi 
Jo 

a\v\ dv\ 



(2.10) 



a±vi dvt 



where, for simplicity, we set a(r, Vx, mi, t) — a\. The 
maximum distance / to the perturber stars has been taken 
equal to the radius of the cluster. The subscript p reminds 
us that the variation of a is due to perturbations. 

2.4 Global evolution equations 

It is now time to focus on the entire cluster. Let's consider 
the subset of stars of given mass m and energy E, at time 
t. These stars are situated at different locations in the 
cluster; but, according to our hypotheses, they all corre- 
sponds to the same value of /, thus of a. And yet, after 
a very short time, the functional form of a has changed. 
This modification, given by (2.10), is different from one 
position in the cluster to another, in general. As a con- 
sequence, a does not take the same value over the entire 
subset of stars anymore. In other words, the perturba- 
tions instantaneously destroy the steady-state regime. 

However, the steady-state re-appears almost immedi- 
ately, because of the rapid circulation of the stars within 
the cluster; the perturbations become equal for all the 
stars of the subset, so that a takes the same value for all 
of them, once again. 

Thus, the actual variation of a, after the equalization, 
is obtained when computing the derivative of (2.9), i.e.: 

da _ d£dU df 
dt ~ dE dt + dt' 



(2.11) 



2 This equation has recently been cited and used by I Kind {i960). 
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and we write that da/dt, given in (2.11), is equal in av- 
erage to (da/dt) , given in (2.10). 

Let g(E,m,t) dE dm be the total number of stars in 
the subset. One has: 



g(E, m, t) dE dm 



Anr 2 dr ■ 4irav 2 dv dm, (2.12) 



where the integral is with respect dr. The maximum value 
r rn taken by r for a given energy E, is set by: 



U(r m ,t)=E, 
and the velocity v is given in (2.3), i.e.: 
v = (2E-2U) 1/2 . 
Equation (2.12) becomes: 



a = 16n / r dr ■ av, 
Jo 

and the number of stars in the subset changes as: 

da d(av) ( dv\ 



m =16lT 



r 2 dr 



dt dv \dt J 



(2.13) 



(2.14) 



(2.15) 



(2.16) 



We write that this integral remains the same whether 
using (2.10) or (2.11). The second term in the integral is 
the same in both cases; therefore we get: 



r drv 



da 

~dt 



d£dU df_ 
dE~dt ~ dt 



0, (2.17) 



which is the equation of evolution we were looking foip. 

By changing a and a\ into / and fi and using (2.14), 
it becomes: 



= J r 2 drl 167rGMn(n)y m 1 dm 1 (2.18) 



d_ 

dE 



mf f f 1 (2E 1 -2U) 1/2 dE 1 
Ju 

(2E-2U) 1 



1/2 fd£dU df 
dE dt dt. 



We can switch the operations J^ m dr and ^ , because 
the quantity between square-brackets becomes zero for 
T = r m . Then, we switch the order of the integrations 
over r, mi and E\, by recalling that / and fi do not 
depend on r. Furthermore, we set: 



{2E-2U) 3/2 r 2 dr = q(E,t) 



(2.19) 



o 



3 This equation has already be en presented, in a slightly different 
form however, by Kuzmin (1957, Equation 2.5). 



([T.N.: The upper limit of the integral is missing in the 
original version]) and thus, 



dq 
dE 
dq 

In 



{2E-2U) 1/2 r 2 dr, 
dU 



(2.20) 



(2E - 2U) 1 ' 2 ^-r 2 dr. (2.21) 



After these transformations, (2.18) takes its final form: 



16tt z G 2 ln(n) 



mf 



f f 1 q' 1 dE 1 +m 1 f'( f fxqxdE 

J — oo J — CO 



/"CO 

+q / fx dEx 

J E 



dt 





d 


dm\ 






~d~E 








r 




-co 




df 




dt 



(2.22) 



(Hereafter, the symbol ' indicates the partial derivative 
with respect to E: §^ — /', and so on.). 

We note that the structure of the cluster is present in 
this equation only via the function g, defined in (2.19). 
This function will play a very important role in the fol- 
lowing. 

We now have the complete system of fundamental 
equations that rule the structure and the evolution of the 
cluster: these are the four equations (2.6), (2.7), (2.19) 
and (2.22), linking together the four quantities /, p, U, 
q. But before actually solving this system, we first have 
to modify it. 

2.5 Equal masses 

In the rest of this paper (except in Chapter VI) , we limit 
ourselves to the study of the case where all the stars have 
the same mass. This is our second arbitrary hypothesis, 
less critical than the first one however. Indeed, whereas 
it is obvious that the masses of the stars are different in 
real clusters, one can imagine a fictitious cluster where 
all masses would be equal. In other words, the system 
we study corresponds to a physically plausible situation, 
while being simpler than the real systems. 
Therefore, the distribution function reads: 



f(E, m, t) = 8{m -mi)- F{E, t), 



(2.23) 



where mi is the mass of each star, and 8 is the Dirac's 
distribution [T.N.: Dirac delta function] . Performing the 
integrations over m in (2.6) and (2.22) is immediate. 

2.6 Normalized variables 

In order to get rid of the numerical constants, we define 
"normalized variables" : 



P 
r 

Q 
dt 



47T7711 D, 

(16tt 2 G mi)- 1 ' 2 R, 
(16tt 2 G mi)" 3 ' 2 Q, 
[16tt 2 G 2 mf ln(n)]- 1 dT. 



(2.24) 
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The four fundamental equations become: 



/>oo 

D = (2E - 2U) 1/2 FdE, 

Ju 



(2.25) 



d 2 U 2 dU 



dR 2 RdR 



1 



= D, 



3 Jo 
_d_ 

dE 



Q 



+Q / F 1 dE 1 
Je 

In addition, we set 



(2E~2Uf /2 R 2 dR, 
F f PlQi d^i + W / FiQa cLEi 

</ — OO J — OO 



F ,dQ_ ,dF 
dT v 9T' 



i^id^i. 



(2.26) 



F 1 Q[dE 1 -F F^ 1) Q' 1 dE 1 . (2.27) 

-oo ./ — OO 



We have 



as 
as 



F 1 Q' l dE l -F' F^Q'idEu 

-oo J — OO 

(2.28) 

i.e., thanks to a partial integration and by taking into 
account that Q(— oo) = 0: 



S' = F F 1 Q[ d.Ei 

./ — OO 



(2.29) 




so that (2.25d) can be re- written: 

()-S" + F> dQ Q' dF 



(2.30) 



2.7 Canonization 



We know (see e.g. iKurthl 119551) that every theoretical 
cluster model includes two dimensional parameters; or, 
in other words: one can apply to the model an homol- 
ogy based on two arbitrary parameters. This can be eas- 
ily checked with the equations: at a given time To, the 
cluster is completely defined by the distribution function 
F(E,Tq). The most general homologous transformation 
would be to multiply F on the one hand, and E on the 
other, with two arbitrary constants. (We assume that 
the masses of the stars are set; thus, it is not possible to 
change them through a homology.) 

It is interesting to use this possibility of homology 
to set any cluster model to a canonical form defined by 
two chosen conditions. This allows one to highlight the 
very structure of the model, in a way that is indepen- 
dent of its size. The two conditions could be, for exam- 
ple: masses and radius equal to defined values. However, 



we will choose others, more convenient (see Chapter IV, 
Equation 4.6). 

In the case of a model in evolution, which we study 
here, this model will always be set to the canonical form 
thanks to a proper homology. Therefore, the two param- 
eters of the homology will be functions of time. This 
has the great advantage to allow us to split two aspects 
of the evolution of the model: (1) the evolution of the 
size, represented by the variation of the parameters of 
the homology; (2) the evolution of the structure itself, di- 
mensionless, represented by the variation of the canonical 
form. 

Therefore, we set: 



E 
F 



(3 E 
7 F, 



(2.31) 



where f3 and 7 are the two fundamental parameters of 
the homology, chosen so that the new model, defined by 
F(E), is canonical. By replacing this in the equations, 
we find out that the transformation formulae of the other 
quantities are: 



(2.32) 



u 


= pv, 


D 


= /? 3/2 7D, 


R 


= p- 1 '* 7- 1 / 2 R, 


Q 


= /? 3 / 4 7 ~ 3/2 Q, 


dT 


= 7 _1 dT, 


S 


= /3 y / 4 7 1/2 S. 



One can easily check that the first three equations of 
the system (2.25) remain the same with the new variables. 
The fourth one changes however, because of derivatives 
with respect to time. From (2.31a) and (2.32e), we get: 



9(E,T) 



8(E,T) 
which allows us to compute: 



f3 2 dT 

7 



d£dQ _ dQdF _ 7 
~dE dT ~ dEdf ~ /? 



dFdQ 
ffEdT 



dQ_dF_ 
<9E dT 



By use of (2.31b) and (2.32d), this becomes: 



dFdQ _ dQdF 
dEdf ~ dEdf 



(2.33) 



(2.34) 



(2.35) 



dFdQ 
dE dT 



dQdF 
dE dT 

2 7dT^ 



Q 



dF ,3 1 d/? 
9E I^dT 



p dQl d7 
dE 7 dT 



Let: 



1 dP 
/?dT 
1 d 7 
7dT 



(2.36) 
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With the new variables, the system (2.25) becomes: 



poo 

D = / (2E-2U) 1/2 FdE, 
Ju 



(2.37) 



d 2 U 2 dU _ 
<9R 2 + R 5R ~ ' 



Q = -J (2E - 2U) 3/2 R 2 dR, 

FiQ[ dEi + F' ( / FiQidEi 

OO \ J —OO 

+Q y Fi dEi^j 

. S» + F'§9-Q<§| + g 6 -5c)QF-cPQ'. 

(The symbol ' indicates here the derivative with respect 
to E, and not E any more.) 

2.8 Secondary equations 

We write here the equations that allow us to compute 
several interesting quantities. 

2. 8. 1 Partial mass and partial kinetic energy 

Let Ve be the subset of stars whose total energy is less 
than E. The total mass of this subset Ve is, when using 
(2.14): 

M= mdm 4irr 2 dr 4ir(2E - 2U) 1/2 .h dE u 
Jo Jo Ju 

(2.38) 

and its total kinetic energy is 

C= mdm 4irr 2 dr 2n(2E - 2U) 3/2 fx dE x . 
Jo Jo Ju 

(2.39) 

Switching the order of the integrations and using 
(2.19) and (2.20), one obtains: 

poo pE 

M = 16tt 2 / mdm / /i^dEi, (2.40) 

JO J-oc 



C = 24tt 2 / mdm 



/ mdm fiqi 

Jo J -oo 



dEi. 



We assume all the masses to be equal and we introduce 
the normalized variables defined as: 

M = (16tt 2 mi y 1/2 G- 3/2 M, (2.41) 
C = (16tt 2 mi y 1/2 G~ 3 / 2 L, 

which results in: 

r-E 



M 

L 



/Hi 
FxQ\ did, 
-oo 
o r-E 

= -z / FiQid^L 

Z J -oo 



(2.42) 



These relations allow us to give a physical meaning to 
the integrals found in the evolution equation (2.25d). We 
also note that if we set E = +oo, the subset Ve matches 
the entire cluster and the equations (2.40) and (2.42) give 
the total mass and the total kinetic energy of the cluster. 

We switch to the canonical variables by setting: 

M = (3 3/4 7~ 1/2 M, (2.43) 
L = /3 7 / 4 7 " 1/2 L- 

2.8.2 Another partial mass 

The mass within the sphere of radius r is 

M r = J 4nr 2 Pl d n = -r 2 — , (2.44) 

which becomes, after transformations as in (2.41a) and 
(2.43a): 

/• R ,9TT 
M R = y DiR^dR^R 2 — (2.45) 

One must be careful not to mix up Mr, the mass within 
the radius R, with M, the mass of the stars whose energy 
is less than E. 

2.8.3 Projected density 

The projected density of the cluster, i.e. the mass ob- 
served per unit surface, is 

p (Vr 2 +z 2 ) dz. (2.46) 

-co 

We switch to the normalized, and then canonical, vari- 
ables with 



PP {r) = m{ /2 G- 1/2 D 



(2.47) 

D P (r) = /3 5 / 4 7 1 / 2 Dp, (2.48) 
and we obtain 

D P = / D (Vr 2 + Z 2 ) dZ. (2.49) 



2.8.4 Projected mass 

The mass enclosed within a circle of radius r in the pro- 
jection of the cluster is: 



Jo 



Mp = I 2nri p P dri, 



(2.50) 



i.e., through canonical variables and transformation as in 
(2.41a) and (2.43a): 



M, 



= — f R D 

2io 



RdR. 



(2.51) 
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2.9 Homologous evolution 

The system of equations (2.37) allows us to compute the 
evolution of the cluster from a given initial state. How- 
ever, a few tests revealed that even with a powerful ma- 
chine (IBM 704), this calculation is extremely long. Fur- 
thermore, the circumstances of the birth of a globular 
cluster (and even more, the initial energy distribution) 
remain unknown at present time. The choice of the func- 
tion describing the initial state of the cluster is, therefore, 
quite arbitrary. 

Consequently, it seems more rational to face the prob- 
lem from the other side, and to first look for what the final 
evolution of the cluster will be. An analogy with gases 
suggests indeed that the stellar systems should tend to- 
ward a defined final state, independent of their initial 
state; and the simil arity of the globu lar clusters as we 
observe them today (jvon Hoernerlll957l ) seems to confirm 
this assumption. This final state should be deduced from 
the system of equations (2.37) alone; it is the equivalent, 
in a way, to the Maxwellian distribution of velocities for 
a gas in equilibrium. 

However, it is well-known that the final state of a clus- 
ter cannot be equilibrium. The stars which gain enough 
energy through perturbations escape forever from the 
cluster, whose total mass therefore decreases. Hence, a 
stationary state stricto sensu, does not exist. 

How then, can we imagine the final evolution of the 
cluster? It is time to go back to the distinction made 
above between the evolution of the size and the evolution 
of the structure. The escape of stars only implies the 
evolution of the size; it does not forbid the existence of a 
model for which only the dimensions would vary, while its 
structure would remain the same; in other words, a model 
where the evolution would only consist in "expansions" or 
"contractions" of the physical quantities. We call such a 
process homologous evolution. 

Therefore, it seems natural to seek whether the final 
evolution of the cluster would be of this kind. This is 
what we will do in the next Chapters. The answer will be 
yes: we will see first (Chapters III to V) that the system 
(2.37) allows one and only one solution in homologous 
evolution, which we will name homologous model; second 
we will show (Chapter VII) that this model represents 
indeed the final state toward which all the clusters tend. 

Writing of the above in mathematical terms is very 
simple. A model in homologous evolution allows an in- 
variant canonical form: F, Q, and so on, are independent 
of the time. Thus, in particular: 



dF 
9T 



= 0, 



dQ 
dT 



= 0, 



(2.52) 



Before focussing on the numerical solution, we still 
have to examine the boundary conditions. This will be 
the topic of the next two Chapters. 

3 BOUNDARY CONDITIONS 
3.1 Effect of the galactic field 

Globular clusters are not isolated systems: they are sub- 
ject to the gravitational field of the galaxy, an effect that 
cannot be neglected. 

Let Uq be the galactic potential. The dimensions of a 
cluster being much smaller than those of the galaxy, this 
potential can be described with sufficient precision by a 
Taylor series truncated to the second order; in a reference 
frame centered on the cluster, moving with it, and with 
the right orientation, the series reads 



U G = U G (0) + - 



1 fd 2 U G 



d 2 U G 2 

y 



dx 2 dy 2 

Furthermore, we almost have: 

d 2 U G d 2 U G d 2 U G 



3 2 U G 
dz 2 



dx 2 



dy 2 



dz 2 



0. 



(3.1) 



(3.2) 



because the density of the galaxy is very small in the 
regions where the globular clusters are found. As a con- 
sequence, at least one of these three main curvatures is 
negative. We assume that the one along x is the most 
negative. 

The potential created by the cluster is, in the external 
regions: 

Ua = (3.3) 
r 

[T.N.: the subscript A stands for amas: cluster.], where 
A4 e is the total mass of the cluster. Thus, the total po- 
tential along the a;-axis is (Figure [1} : 



2 dx 2 



It reaches a maximum for 



x e = (G M e ) 1/3 



U e = U G (0) - ^(G M e ) 2/3 (- 



d 2 u G 

dx 2 

d 2 U G 
dx 2 



-1/3 



1/3 



(3.4) 



(3.5) 



It follows that every star that goes beyond this point 
escapes from the cluster forever. This implies that the 
stars that have a radial, oscillatory, trajectory along the 
x-axis cannot have a total energy greater than U e ; there- 
fore, for these stars: 



and the last equation of the system (2.37) shrinks to 



/ = 



for E>U P 



(3.6) 



^b- ~c) QF'-cFQ'. (2.53) 



This way, the system becomes independent of time. 



But, as we supposed that the distribution function de- 
pends only on E, (3.6) is, in fact, a general property: no 
star can have an energ y higher U e . In U e , the distribu- 
tion function has what IChandrasekhar (1943d) calls an 
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Figure 1: Galactic potential Ug, cluster potential Ua and total 
potential Ug + Ua- 



"absorbing barrier" , and that we call here, in a more il- 
lustrative way, a "leak". We know that, in this case, the 
function / is zero for U e , but its derivative is not; the 
value of t he latter is proportion al to the escape rate of 
the stars (| Chandr asekharll 1 943dl Equations 25 and 33). 

The galactic field has also the effect of modifying the 
potential inside the cluster. However, the equations (3.1) 
and (3.2) show that Ug — Ug(0) is zero in average on a 
sphere centered on the cluster; as a consequence, in our 
approximation of spherical symmetry, the galactic field 
has no effect inside the cluster. Its effect reduces to the 
creation of a "leak" at U e . 

The boundary of the cluster is defined by the equa- 
tion U = U e ; it is a closed surface, with a shape of a 



lemon, elongated along the x-axis (on which it has two 
conic points for x = ±x e ). Therefore, the "radius" of the 
cluster is a function of the angle; but in the case of the 
spherical symmetry approximation, we will limit ourselves 
to considering its mean value. It is obtained from (3.3) 
and (3.4), by replacing the galactic potential Uq with its 
mean value f7c(0), which gives out 



G M e 



U G (0) 
2 
3 



(GM e ) 1/3 



(3.7) 



d 2 u G 

dx 2 



-1/3 



This is this mean value that we will refer to as the radius 
of the cluster, from now on. 

The relation (3.7) links the radius to the mass of the 
cluster. The curvature of the galactic field can be consid- 
ered as constant: we know that the dynamical evolution 
of the galaxy is much slower than that of the clusters. It is 
true that the cluster travels along an orbit in the galaxy; 
but this motion is much faster than the evolution of the 
cluster, and we can use the mean value of the curvature 



of the field along the orbit in (3.7). 
becomes 

r e oc Ml /3 . 



Thus, the relation 



(3.8) 



Isolated cluster 

For future applications, we will also consider the case 
of an isolated cluster, which would experience no external 
gravitational field. Indeed, it seems that several systems 
(galaxies, galaxy clusters) are almost in this situation. In 
this case, at a sufficiently large distance from the center, 
the total potential reduces to Ua, given in (3.3). It does 
not yield a maximum but rather increases with distance. 
Therefore, escapes can occur only if the distribution func- 
tion ranges up to E — 0, corresponding to an infinite ra- 
dius for the cluster. Furthermore, we find out from (2.20) 
and (3.3): 

y- * G * M * (39) 
8V2 {~Ef/ 2 ' 

In order to have a finite total mass (which is given by 
Equation 2.40), the distribution function must decrease 
steeper than (— E) 3 / 2 for E — > 0. It follows that its deriva- 
tive is zero at the leak point E = 0, leading to an escape 
rate of zero. This result has al ready be shown by means of 
another method ( Henon 1960f) : it is confirmed by th e be- 
havior of the artificial clusters o f lvon Hoernen (|l960h from 
which no stars escape, even after a time much longer than 
the relaxation time. Therefore, for an isolated cluster, one 
has: 



M e = est. 



3.2 Mass-radius relation 

The two relations (3.8) and (3.10) can be represented by 
the unique form: 

(3.11) 



(3.10) 



with A = 1/3 for globular clusters (and for any system in 
a non-uniform external field in general) and A = oo for 
isolated clusters. We will use this general form (3.11), so 
that the formulae obtained could be eventually applied to 
the case of isolated clusters; but all the numerical calcu- 
lations, in the present paper, will be done for A = 1/3. 

Using (2.32c), (2.43a) and (2.36), we switch to the 
canonical variables, and the relation (3.11) becomes: 



din (Re) 
dT 



= A 



36 

T 



din (M e ) 
dT 



(3.12) 



If the model is in homologous evolution, R e and M e 
are independent of time, and we get: 



0= (3A + l)6+(2-2A)c. 



(3.13) 



We see that the mass-radius relation translates into a 
definite relation between the two parameters of the ho- 
mology. 
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3.3 Choice of a reference level for the potential 

The potential U is defined up to an additive constant. We 
set this constant by writing: 

U e = 0. (3.14) 

Therefore, the potential is zero at the boundary of 
the cluster, negative inside, positive outside. The total 
energy of a star is always negative, E — corresponding 
to the escape threshold. 

4 CENTRAL CONDITIONS 
4.1 Infinite central density 

A priori, one could expect the density profile of the model 
we seek to have the classical shape of Figure [5] with, in 
particular, a well-defined central density; and our first 
attempts have been set to find this. However, the calcu- 
lations of the evolution from such an initial state revealed 
an unexpected phenomenon: no matter how big the initial 
value of the central density is, it always increases, without 
limits. As an example, Figure [5] shows the spacial density 
computed at two successive times. The horizontal axis is 
normalized to the external radius; this way, we see that 
the increase only concerns the central region of the clus- 
tefl This increase cannot be diminished by means of an 
homologous transformation; indeed, one can measure it 
by computing the ratio between the central density and 
the mean density of the cluster, and this ratio does not 
change in an homology. 

Therefore, we have been led to consider, for the final 
state of the cluster, a model with infinite central density. 
This can seem physically absurd at first sight. However, 
one should not forget that the density considered here 
is a probability density. The cluster is not made of a 
continuous medium, but rather of separated particles. For 
the model to have a physical meaning, the integral of the 
probability density must be finite over any finite volume; 
and we will see it is indeed the case. 

In order to better emphasize this important point, we 
have created, from a table of random numbers, an "arti- 
ficial cluster" , that follows the projected density law 

pp oc -. (4.1) 

r 

This cluster, shown in Figure [31 looks like a real galac- 
tic cluster; the fact that its theoretical density is infinite 
in its center is not visible, and if we would plot the density 
profile from Figure [3l as one would do for real cluster, we 
would likely state that the curve shows a finite maximum 
in the center! 

The functional form of (4.1) is indeed the one we will 
find near the center of the model (see Equation 5.9). Note 
that the number of stars inside a circle of radius r is 

iMichid Jl96lfl has also observed a flux of stars toward the center 
of his model. 



t 1 1 1 r 




Figure 2: Increase of the density in the central region of a cluster 



proportional to r, and thus tends toward zero when the 
circle becomes smaller and smaller. 

In Chapter VIII we will study in detail how this central 
singularity of the density arises. For the time being, we 
only focus on the final state, and thus we assume that the 
singularity exists. 

4.2 Asymptotic expressions near the center 

p being infinite at the center, the function / must become 
infinite too, as shown in (2.6). As a consequence, the first 
two terms of the equation (2.53), which are of the order of 
FQ, become negligible with respect to the first one which 
is of the order of F 2 Q, and we have: 

S" = 0, (4.2) 

i.e., after integration, 

S = ai+a 2 E, (4.3) 

where ot\ and cii are two constants. 

Let's assume first that ot\ and are zero. In this 
case, using (2.27) and (2.28), (4.3) reduces to 



which immediately integrates as 

F = d c- C2E , (4.5) 

where C\ and Ci are two constants, necessarily positive, 
because F must be positive and increasing for E — > — oo. 
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Figure 3: Artificial cluster with infinite central density. 



This way, we find a Maxwellian distribution of the 
energies near the center, independently of the particular 
functional form of Q, i.e. of the structure of the cluster. 
This result is in agreement with what one could expect: 
in the central regions, which are very dense, the perturba- 
tions are very active and almost establish the Maxwellian 
equilibrium^. We are now to set the two conditions that 
define the canonical form of the cluster (see Equation 2.31 
and the following): we assume that E and F have been 
normalized so that the constants C± and C 2 equal unity. 
The distribution function near the center is then: 



F = c" 



(4.6) 



This choice of the canonical conditions has the advantage 
of simplifying the formulae and the calculations a lot. 

From (4.6), using successively (2.37a), (2.37b), 
(2.37c), we easily find the expressions of the various quan- 
tities near the center: 



1/4 



R 



(4.7) 



3 U/2 



Q = K D e 3E / 2 , 

K D = _ 16 -=0.45841). 

9^ 87T)V4 J 



Note that (4.7b) is obtained as the unique solution of the 
differential equation (2.37b) for which D is infinite at the 



5 often called "isothermal equilibrium", which is a misuse of lan- 
guage, because the concept of temperature is meaningless in stellar 
dynamics. 



center of the cluster. We have, near the center: 
Also note that the potential is, near the center: 



U = - In 

2 



2 In (R) 



(4.8) 



(4.9) 



Therefore, it tends logarithmically toward — oo in the cen- 
ter: the cluster yields an infinite "potential well". As a 
consequence, the energy of the stars can take any value 
from E = — oo to E = 0. 
Finally, we have 

/E 
FiQidEi = 2K D e E ^ 2 , (4.10) 
-oo 

/E 
FiQidEi = 3K D e E/2 . 
-oo 

Let's go back to (4.3) in the general case where ct\ 
and a 2 are not zero. For small values of a\ and a 2 , the 
solution of (4.3) can be expanded as: 

F = c- E + aiyi+a 2 Y" 2 , (4.11) 

where Y\ and Y% are function yet to be determined. By 
setting this in (4.3) and using (4.7c), we obtain two dif- 
ferential equations for Y\ and Y 2 : 

e -E/2 

2Y1+F/ = — — , (4.12) 



Y 



(-i) 



Y 



(-i) 



2Y 2 + Y{ 



which can be solved as 

Y x = Cne- E + Ci 2 E c- E 



3^15 

(1 + E) e- E / 2 
3K D 



3 -E/2 



Y 2 = Czie" 



C 22 E 



3K D 
2 



(4.13) 
(5-E) e- E / 2 . 



The Cij are arbitrary constants. The terms Cy are, 
in fact, parasite solutions coming from a bad normaliza- 
tion of F(E); we can get rid of them by modifying this 
normalization. Then, F reads 
2 

(5a 2 -ai-a 2 E)e~ E/2 . (4.14) 



3K D 



The second term of this expansion becomes negligible 
with respect to the first one, for all a\ and a 2 , when E is 
small enough. As a consequence, the expansion (4.14) is 
valid near the center, not only for small a\ and a 2 , but 
for any values they can take. Thus, we have found the 
general solution of (4.3), i.e. the functional form of the 
distribution function, for E — > — ooQ 



6 The expression (4.14) can be considered as an expansion of F 
up to the second order in e E / 2 ; however, it has been obtained using 
the expansion of Q up to the first order, (4.7c). This method is 
valid because S cancels when F takes the form (4.6), for all Q's. 
We easily find out that the expansion of S up to the n-th order 
depends on the expansion of F up to the n-th order, and on the 
expansion of Q up to the (n — l)-th order only. Therefore, it is 
possible to consider one term less in the expansion of Q than in 
those of F. 
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4.3 Flux of matter and flux of energy toward the 
center 

We are going to calculate the flux of matter and the flux of 
energy through the surface E = est, in the six-dimensional 
space of positions-velocities. This calculation will allow 
us to give a physical meaning to the constants at\ and 
a 2 , and to specify the conditions the model must fulfill 
near the center of the cluster. Let's go back a second to 
the study of the local properties and let's consider, in one 
point in the cluster, the subset of stars whose velocity is 
smaller than a given v. The density of this population is 



Pv 



/>oo r-v 

I to dm / Airai vfdvi (4-15) 
Jo Jo 

r-v 

Jo ' 



= mi I AttAi v 1 dvi, 



where a is the local distribution function, defined in (2.9). 
In analogy with (2.23), we set 



a = S(m — mi) A. 
From this, and using (2.10), we find 

= 64tt 3 G 2 ml In (n) 



dpi 
dt 



pV 

7 



A Ai vf dvi 



io 3dv\vJ 

■j 2 I Ai vi dv 
Jo 



(4.16) 
(4.17) 

i 

Ai v\ dvi 



This quantity can be physically interpreted as the flux 
of stars, in velocity-space, through the surface v = est. 

The total mass of the stars in the subset Ve (stars 
whose total energy is less than a given E) is (see Equation 
2.38): 



M 



- r 

JO 



Aup v r dr, 



(4.18) 



which becomes, after differentiation while keeping E con- 
stant: 



dM" 



dt J, 



47T 



dv\ 



n dr. 



d£v\ dpv_ 
dt ) p dv V ^ /B j 

(4.19) 

This is the flux of the stars through the surface E = 
est, in the six-dimensional space. The subscript p of the 
first term can be removed, because the equalization of 
the perturbation occurs between stars of same energy, and 
does not affect the value of M. Furthermore, from (2.14): 



-(2E-2U) 



-1/2 



dU 



(4.20) 



Substituing into (4.19), then changing the variables, 
switching the operations the same way as in Chapter II, 



and using (2.19), (2.20) and (2.21), we obtain: 

,dq 
dt 



i)M : 16^ 2 + 2567T 4 G 2 m 2 In(n) (4.21) 



/E / pE 

Fi q^dEi+F'i / Fi qidEi 
OO \J—QO 



dt 



+q / FidEi 

J E 

We notice that when using (2.40a), we can easily go 
from this relation to the fundamental equation of the evo- 
lution (2.22), which is thus proved again. On the other 
hand, it is not possible to go from (2.22) to (4.21), be- 
cause an integration constant remains undefined; that is 
why we had to go back to the basic equations to establish 
(4.21). 

We switch to normalized variables thanks to (2.24) 
and (2.41a), which gives, using (2.29): 



dM 
~dT 



dT' 



(4.22) 



Let's compute, in the same way, the flux of energy 
through a surface E = est. At a point in the cluster, the 
subset of stars which velocity is less than v has a total 
energy density 



h v =mi 

Jo 



4ttA Ev 2 dv, 



whose variation is 

dh v 
~T)t 



47T 



A dU 
~dt 



v 2 dv. 



(4.23) 



(4.24) 



(4.25) 



We will not give the details of the calculation, which is 
similar to the previous one, but slightly longer; one must 
do several integrations by parts. We set, in analogy with 
(2.41b): 

U = (16^ 2 mi)- 1 / 2 G- 3/2 H, (4.26) 

and we obtain the expression of the flux of total energy 
through the surface E — est: 



The total energy of the subset Ve is 

H = ! 4:irh v r 2 dr. 



dH 



dQ 



= ES' - S + EF r 
dT dT 



£«§£d*. (4.27) 



It is interesting to note that by computing the deriva- 
tive (4.22) and (4.27) with respect to E, we find: 



dH' „ dM' 
= E 



dT 



dT 



(4.28) 



This relation can be proved more directly. Indeed, let 
dVE be the subset of stars whose total energy is between 
E and E + dE; the mass of this subset is: dM = M'dE, 
and its total energy is: dH — H'dE. Thus comes 



H' = EM', 



(4.29) 
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and the relation (4.28). However, here again, it is not 
possible to use (4.22) and (4.28) to demonstrate (4.27), 
because an integration constant would remain undefined. 

(The quantity W, that we called "total energy of the 
subset Ve", has been computed by simply adding the 
total energies E of its members, which is not correct, be- 
cause by doing this, we count twice the mutual potential 
energy of the stars of Ve- However, this effect is negligi- 
ble if Ve only counts a small fraction of the stars of the 
cluster, which is the case for E — > — oo.) 

Near the center, one can neglect the last term of 
(4.22), and the two last terms of (4.27), which are of the 
order of e E / 2 ; we also replace S with its value, given in 
(4.3) and (2.32f), and we get 



dM 
~dT 
dH 



= /3 3 / 4 7 1/2 « 2 



(4.30) 



vV2 



ax. 



Thus, the two constants ot\ and ct2 represent, up to a 
factor, the fluxes of mass and of energy toward the center 
of the cluster. This result is unexpectedly simple. 

Note that we can switch to the canonical variables 
thanks to (2.32e) and (2.43); H has still the same dimen- 
sion as L. By neglecting once more the terms in e B / 2 , we 
obtain: 



dM 
dT 



<*2, 



(4.31) 



-eti. 



We are now going to examine the physical meaning of 
these fluxes. Let's focus first on the flux of mass. (4.30a) 
shows that, near the center, it takes a constant value, 
independent of E\ this means that a flow of matter enters 
or exits the cluster (depending on the sign of a^) through 
the central singularity. The fundamental equations do not 
rule out such a flow; one can even, as we will see, create 
an infinite number of models obeying the equations and 
yielding a non-zero flow of matter in the center. However, 
from a physical point of view, such a flow obviously does 
not make sense. This leads us to write the additional 
condition 

a 2 = 0. (4.32) 

The need of writing a separate condition for the mass 
conservation in the center can be more easily explained: 
at the center of the cluster, the quantities become infinite 
and the fundamental equations are meaningless; thus an 
additional condition is required for this particular point. 
For that matter, we have already noted, just above, that 
it is impossible to go directly from the equation of evolu- 
tion (2.22) to the Equation (4.21) which states the con- 
versation of the mass; one constant is missing, which is 
precisely the value of the central flow. 

However, the models with a non-zero central flow are 
useful: we will see them again in Chapter VII, while 
studying the initial stages of the evolution, during which 



the central density slowly increases until it becomes infi- 
nite. Then, the flow of matter toward the center exists 
and simply corresponds to the slow "filling" of the central 
part of the density profile (see Figure [2]). 

Let's consider now the flux of energy, given by (4.30b). 
This flux is also constant near the center: therefore the 
center of the cluster creates or absorbs energy. One could 
think that, in analogy with what is true for the mass, the 
flow of energy must be zero in the center. But the numer- 
ical computation shows (see the next Chapter) that the 
previously written condition (4.32) completed the defini- 
tion of the solution of the system of equations, which is 
now unique; and this solution corresponds to a value of 
ct\ that is positive and non zero. Thus, we have to ad- 
mit that a flow of energy toward the center exists; more 
precisely: the center of the cluster absorbs some negative 
energy. 

We will come back with more details to this strange 
and very interesting phenomenon in the Chapter V, and 
we will see how it can be physically interpreted. 

4.4 Follow-up on the expansions near the center 

To prepare the numerical computation, it is useful to push 
further the expansions of the various quantities near the 
center. Taking the relation (4.32) into account and setting 



2ai 



K, 



3K D 

the expansion (4.14) of F reduces to 

F = c- E +Xc- E / 2 . 



(4.33) 



(4.34) 



By using successively the fundamental equations 
(2.37a), (2.37b), (2.37c), we find out the two-terms ex- 
pansions: 



D 



R 



Q 



'7T\l/2 




-u 



2V2 K e~ u/2 ) , (4.35) 



3E/2 



e U/2 _ 2L e u 
8V2 



that make (4.7) more precise. From a remark made above, 
we can compute the expansion of F at the third order by 
putting in (2.53) the expansion of Q at the second order 
only, given by (4.35c). The last two terms of (2.53) must 
now be taken into account. We find: 



p(-i) 



- E +Xe- E / 2 + 0. 



(4.36) 



The third term in the expansion of F is zero. 

We stop the expansions here because the calcula- 
tion of the terms of higher order is much more involved: 
the third term of D is not constant but proportional to 
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(— U) 1 / 2 and its coefficient depends on the entire function 
F. Therefore, we will use the expansions (4.35) for D, R, 
Q. Finally, we note the expansions of the two integrals: 



K 



D 



/E 
FiQidEi =(4.37) 
-OO 

3^3 . 



4\/2 



K 2 e 3E / 2 



FiQi dEi 



2e E/2 L_ Rc E_ VS r2 e 3E/2 



5 THE HOMOLOGOUS MODEL 



5.1 Summary of the equations 

The model that we propose to compute is defined by the 
set of equations and conditions (2.37a, b, c, d), (2.53), 
(3.6), (3.13), (4.34), (4.35b) obtained in the previous 
Chapters. Some transformations are yet necessary to set 
the equations in a form that most favors the numerical 
computation: 



Equations: 

D = 
d 2 Z 



p OO 

/ (2E-2U) 1/2 FdE, (5.2) 
Jv 



dU 2 
R 



v du j 



= z 



Q 





i f (2E-2U) 1 / 2 R 3 dU 
3 J -00 

(f f FiQidEi 

I J — 00 

+F'^y FiQidEi-QF 



Hi.-?. 



FiQ'idEij, 

> J 

= (3A + 1)6+ (2 - 2A)c. 
Boundary conditions: 



F(0) 
F(" 1 )(0) 
F 

R 



0, 

0, 



+ ifc- E / 2 forE^- 



1/4 



,U/2 



7T 



CO, 



for U 



(5.3) 



The unknown functions are: F(E), Q(E), D(U), 
Z(U), R(U), defined from —00 to 0; the unknown con- 
stants are: b, c, K. A equals 1/3. 



1. In the differential equation (2.37b), we will consider 
U as the independent variable. Furthermore, R 
varies rapidly near the boundary, and it is useful to 
switch to a new variable Z defined as: 



R = 



(5.1) 



2. The integral (2.37c) is transformed by means of an 
integration by parts. 



3. The equation (2.53) is replaced with its integrated 
form with respect to E; we have seen in Chapter IV 
that the integration constant is zero. 



All the equations, modified this way, and the bound- 
ary conditions are gathered below: 



5.2 Numerical solving method 

The form of the system does not allow the solution to 
be computed directly; one has to proceed by trial and 
error. Experience leads us to adopt the following iterative 
method: 

1. Choose a temporary form for the function F, that 
fulfills (5.3a) and (5.3c), with a temporary value of 

2. compute D using (5.2a); 

3. compute Z by integrating (5.2b) from the center 
to the boundary; the initial conditions are given in 
(5.3d); 

4. compute Q by means of (5.2c) and (5. 2d); 

5. choose temporary values for c and K; compute b us- 
ing (5.2f); 

6. integrate (5.2e), which is equivalent to a differential 
system of the fourth order, from the center to the 
boundary; the initial conditions are given in (5.3c), 
(4.36a) and (4.37); 
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7. in general, this integration gives the final values of 
F(0) and F(" 1 )(0) that are different from zero; go 
back to point [5.], modify c or K and begin the inte- 
gration again; grope around this way for the values 
of c and K until F(0) and F^ 1 )(0) vanish; 

8. go back to point [2.] with the new function F. 

This way, we obtain a series of approximations for F; 
we stop when two successive approximations are equal, 
up to the desired precision. In practice, the convergence 
is quite fast: the errors are divided by about 5 at each 
iteration. 

The computation has been done thanks to the IBM 
650 device of the Observatoire de Meudon; 8 hours of 
computation are required to get the solution with a pre- 
cision of 1/1000. 



5.3 Results: structure 



R e = 4.703, (5.5) 
M e = 1.996, 
L e = 1.423. 



The four fundamental functions F, D, R, Q are given in 
Table Q] and plotted in Figure |U The table covers the 
range (—5,0); below E = — 5 (or U = —5), the func- 
tions are represented with a sufficient precision by the 
expressions (4.34) and (4.35). The values found for the 
constants are 

K = -0.9499, (5.4) 
c = +0.4078, 

b = -\c = -0.2719. 

Furthermore, the interesting physical quantities take 
the following values: 

external radius: 
total mass: 
total kinetic energy: 



The spacial density D is plotted in Figure [5] as a func- 
tion of the radius. It increases very rapidly toward the 
center, as expected from (4.8). The structure of the model 
is perhaps better rendered in Figure H2 which shows the 
mass fraction Mr enclosed in a sphere of radius R (see 
Equation 2.45). Near the center, Mr is proportional to 
R. We note that half of the total mass is enclosed inside 
the radius R = 0.6800, i.e. only about 1/7 of the external 
radius. 

The potential U vanishes at the boundary of the clus- 
ter, according to our conventions; outside the cluster, its 
form is obtained by integrating (2.45): 

U = Me (R~~R") f°rR>R e . (5-6) 
In particular, for R — > oo, the potential tends toward: 
Uoo = ^ = 0.4243. (5.7) 

tie 

(This is the potential created by the cluster only.) 




Figure 4: Homologous model: the four fundamental functions 
F(E), the distribution function; D(U), the spacial density; R(U), 
the distance to the center; Q(E), see (2.19). 




Figure 5: Homologous model: spacial density and projected den- 
sity as a function of the radius. 
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Figure 6: Homologous model: mass enclosed in a sphere of radius 
R. 



The projected density, computed from (2.49) is given 
in Table [TJ column 6, and plotted in Figure [5] as a func- 
tion of the distance to the center (see also Figure [T5] and 
Figure ITBf . Near the center, from (4.35a) et (4.35b), the 
spatial density is expressed as a function of the radius as: 



D 



_2_ (87T) 1 / 4 K 

R2 + R 



(5.8) 



We find 



D, 



2tt 
Ft 



2(8tt) 1/4 K ln(R) + K P 



for R -> 0, 
(5.9) 

where Kp is a constant which depends on the entire func- 
tion D. By linking the formula (5.9) with the values of 
Table [D one finds K P ~ -3.82. 

The projected mass (mass enclosed, in projection, in- 
side a circle of radius R), computed from the formula 
(2.51), is given in Table [TJ column 7. In particular, we 
find out the value of the median radius Ro of the clus- 
ter, defined as the radius of the circle which contains, in 
projection, half of the total mass. This quantity has the 
advantage of being easily measured for real clusters, while 
the external radius is, on the contrary, almost impossible 
to observe. We find 



Rn 



0.4997. 



(5.10) 



We note that this median radius Ro is about 10 times 
smaller than the external radius R e . 

5.4 Evolution 

When integrating (2.36), we get 



A> e' 
7o e' 



bT 



(5.11) 




Figure 7: Homologous model: evolution of the parameter T, as a 
function of time. 



The "time" T is defined by the differential equation 
(2.32e); it is not proportional to the physical time T. It is, 
somehow, the "proper time" of the cluster; its variation 
is measured with a scale which is proper to the cluster, 
and which always varies according to the evolution. To 
avoid any confusion, it is preferable to consider T as a 
simple parameter which measures the level of evolution 
of the cluster, as it is in (5.11). 

The relation between T and the physical time T, 
found from (2.32e) and (5.11b), is 



1 



ln(l - 70 c T), 



(5.12) 



(by setting the origin of time at T — for T = 0). This 
relation is plotted in Figure [TJ We see that the evolution 
does not last forever, but rather ends abruptly at a time 
7\, given by 

T\ = — . (5.13) 

7o c 

For T = Tx, the mass of the cluster becomes zero, as 
we will see below; it is therefore the time when the cluster 
disappears, after the escape of the last stars. 

However, the curve extends forever in the past: the 
age of the cluster can be anything, it is not possible to 
give it an upper limit. 

When introducing (5.12) in (5.11), and taking (5.4c) 
into account, we obtain the variation laws of the two pa- 
rameters of the homology: 



P = A) i 



7 



7o 1 - 



T 
T[ 

T 



2/3 



(5.14) 



cT 



We derive, from (2.32), (2.43), (2.48), the variations 
of the various physical quantities as functions of time. In 
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particular, for the total mass, we get 



M K = M,, 



T 

1 



(5.15) 



which shows that the mass decreases linearly with time. 
The absolute escape rate is thus constant. 

The radius is proportional to (1 — T/Ti) 1 / 3 , thus de- 
creases quite slowly. The density is constant, and as 
a consequence, the orbital period of the stars within 
the cluster is also constant. The velocities decrease as 
(1 - T/Ti) 1 / 3 . The total energy of the cluster decreases 
(in absolute value) as (1 — T/Ti) 5 / 3 , thus faster than the 
mass. 

5.5 Accumulation of negative energy in the cen- 
ter 

As we have seen in the previous Chapter, a non-zero value 
of K implies the existence of a flow of energy toward the 
center. This energy cannot vanish; we must admit that it 
constitutes an energy reservoir. 

Thus, the total energy of the cluster is always made 
of two parts: 

• a "point" energy Hi, accumulated in the center of 
the cluster; 

• a "diffuse" energy H2, distributed within the entire 
cluster. 

Let's look at how these two fractions of the energy vary 
with time. The variation of Hi equals the flux of energy 
toward the center, i.e., from (4.30b) and (4.33): 



dHi 



£7/4 



1/2 



-K D K. 



dT '' ' 2 
The diffuse energy H2 is, from the virial theorem: 



H> 



-L K 



(5.16) 



(5.17) 



where L e is the total kinetic energy of the cluster; this 
relation remains valid as soon as the potential is set to 
zero at infinity (instead of the convention U e = adopted 
up to now). Therefore, from (2.43b) and (5.2f), we have: 



dH 2 
dT 



dL e 



= -7 



dT 

7 L e 



(5.18) 



1/2 



4-2A 
3A+ 1 



We can finally compute the total variation of energy 
Hi + H2 of the cluster. This variation is only due to the 
fact that the stars escape taking some energy away with 
them. The difference of gravitational potential between 
the boundary of the cluster and infinity is Uoo given by 
(5.7). Therefore, each star that escapes takes away a mass 
mi and a negative energy, equal to —miUao. The escape 
rate is, from (2.43a) and (5.2f): 



dM e 
dT 



7- 



dM e 
dT 



.3/4 -.1/2 



(5.19) 



-^ /4 7 



3A + 1 



M R c, 



and thus, 

d{Hi + H 2 ) 
dT 



£7/4 y/2 



2 M 2 e 
3A + 1 Rl 



c. 



(5.20) 



The comparison of (5.16), (5.18), (5.20) shows that we 
must have: 



3 4- 2A 

2 M+ 3ATI LeC= 3A + 1 R 



2 Ml , 

c. (5.21) 



which expresses the conservation of energy. This relation 
between the parameters of the model provides a useful 
verification of the calculations. When using the numeri- 
cal values (5.4) and (5.5), and when omitting the factor 

£7/4 7 l/ 2) we get . 



d(Hi 



dHi 

~dT 

dH 2 

dT 
-H 2 ) 



dT 



-0.6464, 
-0.9672, 
-0.3455. 



(5.22) 



The conservation of energy is quite well verified; the 
discrepancy that remains comes from the errors in the 
computation (mainly the error made in R e ). 

Thus, in the course of the evolution, the diffuse en- 
ergy H 2 , always negative according to (5.17), decreases 
in absolute value; the numerical values (5.22) show that 
about one third of the negative energy is taken away by 
the stars that escape, while the other two thirds go in the 
center. 

It remains to explain the mechanism of this accumu- 
lation of energy in the center. It is not associated with an 
accumulation of matter; we have supposed that the flux 
of matter toward the center of the cluster is zero (and 
this is indeed necessary, because a central condensation 
of mass would create an additional potential that would 
modify the structure of the cluster; for example the po- 
tential would vary as 1/r and not In (r) anymore, near 
the center). Thus, the negative energy Hi must be stored 
without any increase of the number of stars. Apparently, 
there is only one process that allows this: the formation of 
tight binary or multiple stars in the center of the cluster. 

The direct observation cannot confirm the existence of 
this phenomenon: the images of the stars are sorely sepa- 
rated [T.N.: resolved] in the central region of the globular 
clusters, and it would be impossible to discover there the 
presence of a particularly c ompa ct group of stars. On the 
other hand, von Hoerner (Il960h has computed numeri- 
cally the evolution of artificial clusters; by this means, it 
is possible to observe the mechanism of the evolution, in 
as much detail as desired, von Hoerner has indeed noted 
the frequent formation of binaries in the center of the 
cluster; in one case, a compact group of 4 stars appeared. 
This seems an excellent confirmation of the process that 
we have been led to admit. 

Note again that this central accumulation of energy 
does not affect the structure of the cluster, except for the 
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few stars that support it; thus, it can develop indepen- 
dently of the global evolution of the cluster. In fact, we 
have seen that the absolute value of the central energy 
increases, while that of the diffuse energy decreases. At 
the end of the evolution, about two thirds of the initial 
negative energy is in the central condensation, and one 
third only has left the cluster, carried by the stars. It 
would be very interesting, although not possible here, to 
study in more detail the process of accumulation and to 
answer in particular the following questions: how do mul- 
tiple stars form? How many are there? What happens to 
them after the cluster has disappeared? 



6 STARS WITH DIFFERENT MASSES 

The case of a cluster with an arbitraty mass distribution 
seems much more involved than the simple case of equal 
masses that we have considered so far; thus we shall not 
seek to treat it in the general case. But we arc going to 
see that some simplified cases allow easy calculations, and 
provide, not a complete and rigorous solution, but at least 
several indications about the effect of the dispersion of the 
masses. We will first assume that we add a small number 
of stars of different masses to a cluster that contains equal 
mass stars, and will study the behavior of this secondary 
population. Later, this study will help us to obtain an 
approximate solution to the general case. 



6.1 Simplified model made of two populations 

We suppose that the cluster is made of the mix of a main 
population 1 of stars of mass mi, and of a secondary 
population 2, numerically negligible with respect to the 
first one, of stars of mass m^- The distribution function 
is not (2.23) anymore but rather: 

f(E, m, t) = 5{m - m\)Fi(E, t) + S(m - m 2 )F 2 (E, t), 

(6.1) 

with: 



F 2 < Fi. 



(6.2) 



The population 1 is almost not affected by the pres- 
ence of the population 2; it alone determines the structure 
and the evolution of the cluster; as a consequence, all the 
results obtained in the previous Chapters still work for it. 
We assume that the population 1 has reached the final 
state represented by the homologous model of the Chap- 
ter V. 

Furthermore, by setting m = m 2 in (2.22), one gets 
a new equation, that describes the evolution of F 2 . One 
can neglect the perturbations between the stars of the 
population 2, and the equation becomes, after applying 



the transformation of (2.24): 
d r 



° = 8E 
+F 2 



[iF 2 \ FiQidBi 

J — OO 

E 



(6.3) 



/ F 1 Q 1 dE 1 +Q F 1 dE 1 

J -oo J E 



-oo 

dT 



where we set: 



m 2 

= M- 

mi 



(6.4) 



This equation has the dimension of F 2 ; therefore, we 
can apply an homologous transformation to F 2 , indepen- 
dently of those of F\. Thus, we set 



F 2 — 72 F2 



(6.5) 



Therefore the homology now depends on the three param- 
eters (3, 7, 72. We set 



1 d7 2 
72 dT 



= c 2 , 



(6.6) 



and when continuing the calculation, as in the Chapter II, 
we obtain the equation 







Lf 2 f FiQidEi 

J J — OO 



(6.7) 



rOQ 

F1Q1 dEi + Q / F 1 dE 1 
, Je 



- (V|c)F,Q 



C2 



F 2 QidEi }. 



In this equation, Fi, Q or Qi, b, c are the func- 
tions and constants of the homologous model, given in 
the Chapter V; the function F 2 and the constant c 2 are 
unknown. 

The boundary conditions are only: 



F 2 (0) =0, 



(6.8) 



which tells us that the stars of the population 2 escape 
when they reach the boundary of the cluster. 

Near the center, on the other hand, we can neglect the 
last terms of (6.7) in a first order approximation, and use 
the asymptotical expressions (4.6) and (4.7c) for Fi and 
Q; the equation becomes: 



F' 2 = 0. 



HF 2 

F 2 = C c^ E , 



(6.9) 



and thus 

(6.10) 

where C is an arbitrary constant. The equation (6.10) is 
as expected: indeed, it shows that the Maxwellian equilib- 
rium is realised in the central region of the cluster between 
the two populations. 



20 M. Henon, translated by F. Renaud 



6.1.1 The n > 3/2 case 

This functional form for F2 has a strange consequence. 
The mass of the subset of stars of the population 2 whose 
energy is below E is, from (2.40a): 

M 2 = 16tt 2 to 2 / F 2 q' 1 dE 1 , (6.11) 

J — 00 

i.e., when switching to normalized variables, and then to 
the canonical variables, 



M? 



FaQidEi. 



Q' is given near the center by (4.7c), i.e., 



Q' = ±K D c 3E / 2 . 



5.12) 



(6.13) 



When using this form and (6.10), we immediately see 
that the integral (6.12) diverges if 



(6.14) 



In other words, in this case, the mass of the popula- 
tion 2 is infinite, because of a too rapid increase of the 
distribution function (and of the density) toward the cen- 
ter. 

This anomaly is easy to explain. We have assumed 
in (6.2) that F 2 is always negligible with respect to Fi; 
however, when comparing the asymptotical forms (4.6) 
and (6.10) of these two functions, we can see that when 
/j, is greater than unity, F 2 increases faster than F\ for 
E — > — 00, and that, a value of E below which F 2 becomes 
larger than F\ always exists, however small the constant 
C. This critical value is (from Equations 2.31b and 6.5) 



In (C 72/7) 



(6.15) 



Using the formula (6.13) for Q' is only meaningful 
when E ^> E c . On the contrary, for E <C E c , Fi becomes 
negligible with respect to F 2 . Then, we can obtain the 
new forms of D, R, Q by using (6.10) and the fundamen- 
tal equations. In particular, we get 



Q' 



^V /4 Ik d c 3 ^ e / 2 for E « E c 



7 



(6.16) 



[T.N.: The c subscript is missing in the original version.] 
and with this correct form, the integral (6.12) does not 
diverge anymore. The quantity F 2 Q' yields a maximum 
in the vicinity of E = E c and exponentially decreases in 
both sides. We can compute the mass M 2 in an approx- 
imate way by assuming that Q' is given by (6.16) when 
E < E c and by (6.13) when E > E c ; we find 



M 2 = 3i^(V 5/4 



+ - 



(i \ 7 / C 72 



2[i - 3 / 7 2 V 7 



1/(2^-2) 



(6.17) 



for E > E c , 



M 2 does not depend on E anymore; indeed, because of 
the form of F 2 Q', almost all the stars of the population 2 
have an energy of the order of magnitude of E c , therefore 
are part of M 2 for E » E c . 

We focus here on the extreme case where the pop- 
ulation 2 is negligible with respect to the population 1; 
therefore C must be very small compared to unity. For 
C — > 0, (6.15) shows that E c — > —00. This way, we obtain 
the following result: the stars whose mass is greater than 
3/2 times the mean mass are almost all gathered near the 
center of the cluster. Obviously, this conclusion is related 
to the simplified mass distribution that we have adopted, 
and should not be extended without further study to the 
case of an arbitrary mass distribution. 

We also notice that M 2 is not proportional to C: the 
power of C in (6.17) is less than unity. As a consequence, 
for C — > 0, the last term of the equation (6.7) decreases 
slower than the others, which are proportional to C. In 
the limit, the equation becomes 



-b-c 2 = 0. 



(6.18) 



Let's divide all the terms of (6.7) by C; the last term 
is then, for C — > 0, an indeterminate form x 00. Let 
p be its value; p is a new constant yet to be determined, 
which replaces c 2 . 

6.1.2 Expansion near the center 

Substituting in the equation for the expansions (4.35c), 
(4.36a), (4.37), we get, after some calculations, the ex- 
pansion of F 2 near the center. For < 3/2: 



F 2 = Ce 



(6.19) 



MM-j jj R 2 , c 2 - fic + 1)6 
2 3-2^ 



For n > 3/2: 



F 2 = Cc^fll/ife^ 



+ 



V 

mO - 1) 



(6.20) 



K 2 + lc 



1 



2 2 2 

Finally, for the special case fi = 3/2: 



(6.21) 



-K l + -c-b- 



P 



3K D C 



6.1.3 Method of solution 



Equation (6.7) is integrated from the center toward the 
boundary. The initial conditions are given by one of the 
expansions (6.19) to (6.21). The value of the factor C does 
not matter here; in practice, one sets C = 1. Proceeding 
by trial and error, we find the value of the parameter c 2 or 
p for which the condition (6.8) is fulfilled at the boundary. 
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Figure 8: Stars with different masses: distribution functions. 



6.2 Results: structure 

The calculation has been done for several values of the 
relative mass [i. The value of the parameter is given in 
Table [21 column 2 or 3. Figure [3] plots the ratio of the 
distribution function F2 to its asymptotical form C e _AlE . 
This ratio tends toward unity when E goes to —00; toward 
the boundary, it becomes smaller and smaller, which ex- 
presses the discrepancy between the real distribution and 
a Maxwellian distribution^. 

'Spitzcr Sz Harm ( 1958) have computed these functions in a sim- 
pler case: the structure of the cluster was supposed to remain con- 
stant and was represented by a constant potential inside the cluster, 
zero outside; furthermore, the distribution function F\ of the main 
population was supposed to be Maxwellian. The curves obtained by 
these authors have some similarities with ours but are arranged in 
reverse order! This peculiarity is probably linked to the too rough 
approximation made to the potential. 



-2 1 i 1 ^ — VAXi -I 

1 2 R 3 * 

Figure 9: Stars with different masses: projected density times 
radius, as a function of the radius. 



Figure l^and Table [3] give, for several values of /x, the 
product DpR of the projected density and the radius, as a 
function of the radius; we will see in the Chapter VIII that 
this function is the one that best allows for a comparison 
with the observations. The factors have been adjusted so 
that all the curves cross at the same point: DpR = 1 
for R = 1. The values of C are given in the last row of 
Table El 

For ji = 1, we naturally retrieve the distribution of the 
homologous model; near the center, the quantity D p R 
tends toward a constant (see Equation 5.9). For < 1, it 
tends toward zero; for /z > 1, it tends to infinity. Toward 
the boundary, the curves tend to become parallel. 

6.3 Escape rate 

From the formula (6.11) we find out that the total mass 
M. 2e of the population 2 is proportional to 

/? 3/4 7" 3/2 72- 

As a consequence, the relative escape rate of the stars 
of the population 2 is given by: 

= -J-<^L = -h + 3 -c-c 2 . (6.22) 
M 2e dT 4 2 y ' 

8 is given in Table [2] column 4, and plotted in Figure [TU] 
as a function of /i. From (6.18), we have 

6 = for fj, > 5 (6.23) 

Therefore, in our simplified model, the stars more massive 
than 3/2 time the mean mass do not escape. 



Table 2: 

C2 V 





0.2 
0.4 
0.6 
0.8 
1 

1.2 
1.4 
1.5 
1.6 
1.8 
2 

2.5 

3 

1 



-0.6240 
-0.3998 
-0.1838 
+0.0233 
+0.2208 
+0.4078 
+0.5829 
+0.7430 



0.4796 
0.4670 
0.4307 
0.3859 
0.2685 
0.1726 
0.0685 



1.4396 

1.2154 

0.9994 

0.7923 

0.5948 

0.4078 

0.2327 

0.0726 
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Table 3: log(D P R) 
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log (C) 


1.737 


1.769 


1.797 


1.821 


1.841 


1.857 


1.870 


1.877 


1.858 


1.827 


1.805 


1.769 


1.741 
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Chancirasekh3r,19« 

... Spitzer et Harm, 1858 
— Present modele 




Figure 10: Escape rate as a function of mass, 
reads Present modele: this model. 1 



[T.N.: The legend 



Between fi = and /i = 3/2, the escape rate decreases 
regularly, as one would expect. The very low-mass stars 
(/i ~ 0) escape about 3.5 times faster than the mean-mass 
stars. 

The escape rates computed by IChandrasekharl ( 19421) 
and by Spit zer fc Harm! (|l958l ). normalized so that they 
all take the same value for /x = 1, are also plotted in Fig- 
ure[lJJ We can see that the three curves are very different, 
especially for the low masses. This shows strikingly the 
lack of rigor of the cluster theory, in its present state. 

It is quite obvious, however, that the less massive stars 
must escape faster. The fact that the curve from Chan- 



drasekhar yields a maximum and then decreases for the 
low-mass stars is due to excessive simplifications. Indeed, 
in his computation he assumes that the stars must first 
adopt a Maxwellian distribution (and thus, high velocities 
for the low-mass stars), before being able to escape. In 
reality, as we have seen, the stars escape almost as soon 

as they have reach t heir escape velocity 

Previous works ( van den Bergb 1957 : Takase 1960) 
used the escape rate provided by Chandrasekhar to calcu- 
late the initial mass function of open clusters from their 
present-day mass function; it would be advisable to do 
these calculations again with a more exact escape rate. 
In particular, the abnormal lack of low-mass stars, found 
by Takase in the initial distribution of the Pleiades, is 
likely to disappear. 

6.4 Approximation in the general case 

We are going to show that the results obtained above in 
the case of the very special mass distribution (6.1) allow 
us to solve in an approximate way the general case of any 
mass distribution. 

In general, the distribution function takes the form 
(2.4): f(E,m,t). Let n m dm be the number of stars of 
the cluster whose mass is between m and m + dm; we 
easily derive (see Equations 2.38 and 2.40) 



= 16tT 



fq'dE. 



(6.24) 



n rn depends on m and also on t, because of the escape. 
The total number of stars and the total mass are given 

by 



n m dm, 



(6.25) 



M e = 



i dm. 
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The mean mass and the quadratic mean mass are de- 
fined by means of the classical formulae: 



to 



1 f°° M e 

= — n m dm— , (6.26) 

n J n 



n 



n m to 2 dm. 



This done, let's consider the fundamental equations (2.6) 
and (2.22). We can switch the order of integration with 
respect to to and E. Then, we notice that to can be 
removed if wc introduce the two functions: 



/to dm 



Ji(E,t), 



(6.27) 



f 

Jo 



fm 2 dm = J 2 {E,t). 



Ji and J 2 are "mean distribution functions" , obtained by 
weighting in two ways the distribution functions corre- 
sponding to the different values of the mass. To come 
back to equations already solved, we have to make an 
approximation: wc suppose that theses two mean distri- 
bution functions are similar up to a factor. In this case, 
the equations (6.24) and (6.26) show that J\ and J 2 must 
be proportional to to and m 2 respectively, and thus we 
set 



Ji(E,t) 
ME,t) 



m F(E,t), 

to 2 F(E, t). 



(6.28) 



When using these expressions in the fundamental 
equations, they become 



= 4ttto / (2E -2U) 1/2 FdE 
Ju 



167T 2 G 2 In(n)A 



(6.29) 
(6.30) 



mmf 
+rr?f 



F iq [ dE 1 

E 

F lQl dE x 



+qJ E F^E^j 



■' dt q dt 



The other two fundamental equations, (2.7) and 
(2.19), do not involve the mass, and thus are not modi- 
fied. We now define the normalized variables by means of 
transformations, slightly different from (2.24): 



p = An to D 

r = (16-7T 2 G to)~ 1/2 R 

q = (16tt 2 G to)" 3 / 2 Q 

dt = [16tt 2 G 2 Trf ln(n)]" 1 



(6.31) 



dr. 



We easily check that, this way, we retrieve the equa- 
tions (2.25) which have been obtained for the case of equal 



mass stars. (The equation 6.30 must be multiplied by 
to dm and integrated.) As a consequence, the structure 
of the cluster, given by the functions D, i?, Q, is that of 
the homologous model, and the mean distribution func- 
tion F, introduced in (6.28), also matches that of the 
homologous model. 

Furthermore, when setting 



m 

= m = \i 

m z 



(6.32) 



and replacing / with F 2 , we find that the equation (6.30) 
transforms into (6.3), i.e. the equation obtained at the 
beginning of this Chapter for the simplified model. It 
follows that the detailed distribution function / is repre- 
sented, for the different masses, by the solutions F 2 of the 
simplified model taking (6.32) into account. 

From the previous relation, we note that /x is the rel- 
ative mass computed by considering the mass unit to be, 
not the mean mass to, but rather 



too 



TO 



TO 



(6.33) 



which is different, in practice, by a factor greater than 2 
(see Equation 8.9, below). 

Finally, we note that the transformation equations 
(2.41) and (2.47) must be replaced with: 

M = (16tt 2 to)- 1 / 2 G~ 3/2 M (6.34) 

C = (167T 2 to)- 1 / 2 G-3/ 2 L 
pp = m^G-^Dp. 

7 APPROACH OF THE HOMOLOGOUS 
MODEL 

In this Chapter, we are back to the hypothesis of equal 
masses, and we are going to try to extend the results 
of the homologous model in another direction: we will 
study the evolution of a cluster whose shape is close to 
those of the homologous model, but not identical. We 
will first consider the case of a cluster that differs from 
the homologous model because its central density is finite; 
then the case of a cluster with infinite central density but 
with differences in the global structure. In both cases, 
the calculation will only be approximate. Finally, the 
results will be combined to draw a general picture of the 
evolution of the cluster. 

7.1 Formation of the central singularity 

Let's consider a cluster that matches the homologous 
model everywhere, except in a small central region, where 
it differs so that its central density is finite. The cen- 
tral potential is then also finite; let's call it Uo- Near 
the center, the radius R is proportional to (U — Uo) 1 / 2 ; 
From (5. 2d), we easily derive that Q is proportional to 
(E — Uo) 3 . If we assume that the distribution function 
F also remains finite in the center, we find that FQ' is 
proportional to (E — U ) 2 near the center. It will be 
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Figure 11: Homologous model (dashed line) and simplified model, 
with finite central density (solid line). 



useful, for the next calculations, to have a formula for 
FQ' as simple as possible; therefore, we define a simpli- 
fied model by means of the following conditions: below a 
given value E = Ei, FQ' is proportional to (E — Uo) 2 ; 
above this value, it is as the homologous model. In addi- 
tion, the function and its derivative must be continuous 
for E = Ei (Figure [TT]). These conditions translate into: 



E x = U + 4 
_3_ 
32 



(7.1) 

K D c El/2 (E - U ) 2 for U < E < Ei 



FQ' 



for E > Ei 



We now suppose that the formula (7.1) is valid not 
only initially, but always, although the central potential 
Uo is a function of the time. Thus the evolution of the 
cluster consists in two phenomena that overlap: the nor- 
mal homologous evolution, and the variation of the struc- 
ture of a small central region. We will see that the latter 
can quite easily be computed when considering the fluxes 
of mass and energy toward the center of the cluster. 

The partial mass of the stars whose energy is less than 
Ei is, for the model considered here (see Equation 2.42) 

M = [ 1 FQ' dE = 2Kb c El ^ 2 . (7.2) 

For the homologous model, this mass would be, from 
(4.10b) 

M = 3K D e El/2 . (7.3) 

Furthermore, the mass of the stars whose energy is 
greater than Ei is the same for both models. The differ- 



ence of mass between the present model and the homolo- 
gous one is therefore 



AM = M - M = -K D e El/2 . 



(7.4) 



In the same way, we compute the total energy of the stars 
whose [T.N.: individual] energy is less than Ei: 



H 

H 



EFQ' dE = (6 + 2U )K D c El / 2 (7.5) 



u„ 



(6 + 3U )K D e El / 2 , 



thus, 



AH = H-H = -K D U e El / 2 . (7.6) 

The definition (7.1) of the model is partly arbitrary. 
One can redo the calculations above for other definitions 
(for example, when supposing that Q' rather than FQ' is 
represented by a parabola near the center); we note that 
the results are always 



AM 
AH 



fci e u °/ 2 

(fciUo + fc 2 )e u «/ 2 , 



(7.7) 



where k% and k 2 are two numerical constants, the values 
of which are slightly different from one model to the other. 
In the present case, we have 



ki 

k 2 



-K n e 2 



-3.387 



(7.8) 



0. 



By computing the derivative of (7.7) with respect to 
time, we obtain the fluxes of mass and energy toward the 
center: 



dM = h Uo/2 dUo 
dT 2 6 dT 



(7.9) 



dU 



fjjjS = \{k 1 U Q + 2k 1+ k 2 )e^n 



These fluxes are related, through the equations (4.31), 
to the factors ol\ and a 2 that appear in the expansion 
(4.14) of F. It is more handy to write this expansion as 



F = e~ E + (K + K 2 E) c~ E/2 , 



by setting 



2>K D 



(5q 2 - otx) 



3K D 

Then, we obtain the relations 



«2 



K 
K 2 . 



(7.10) 



(7.11) 



K 
K 2 



fciUo + 7fci + fc2 Un /2 dU 



3K D 



3K D 

a Uo/2 



dT 



(7.12) 



dUo 
dT 



Furthermore, we have seen in Chapter V that the 
boundary conditions require two relations for the param- 
eters of the model. In the case of the homologous model, 
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K 2 is zero, and the other two parameters K and c are de- 
termined uniquely. Here, there are three parameters: K, 
K 2 and c. Therefore, the boundary conditions require, 
after canceling of c, a relation between K and K 2 . As the 
models considered here are close to the homologous one, 
K 2 is slightly different from zero, and we can assume that 
the relation is linear, i.e. 



K = K Q + d Ko 



(7.13) 



where Kq is the value taken by K in the homologous 
model, and d is a constant. 

To find this constant, two models have been computed 
with a non-zero K 2 . The procedure is as in Chapter V, 
but the initial condition (5.3c) has to be replaced with 
(7.10); the other initial expansions have, of course, to 
be modified accordingly. While doing this, we assume 
that the first order expansions (4.7) are still valid; indeed 
these new models differ significantly from the homologous 
model only in the small central region, i.e. for very small 
values of E or U. Near the boundary, the modification 
only affects the second order terms in the expansions, as 
visible for example in (7.10). 

The results about the relation between K and K 2 are: 



K 2 


K 





-0.938464 


-0.01 


-0.931757 


-0.1 


-0.87028 



(K 2 = corresponds to the homologous model; the 
value obtained here is slightly different from that of Chap- 
ter V, Equation 5.4a, because of a larger integration step). 
We verify that these figures match well a linear relation 
of the type (7.13) with 

d= -0.6707. (7.14) 

When combining the relations (7.12) and (7.13), we 



get 



dU _ 3K D K Q c- u °/ 2 
dT ~ hi U + 7 + d + k 2 /ki 
0.3816 e~ u °/ 2 
U + 6.329 



(7.15) 



This differential equation allows to compute the vari- 
ation of the central potential Uo as a function of time, 
and thus to solve our problem. Its explicit solution is: 



T T 2 = 5.241(U + 4.329) e Uo/2 , 



(7.16) 



where T2 is a constant. This relation is plotted in Fig- 
ure [121 We first note that the central potential Uo de- 
creases with time, which corresponds to an increase of 
the central density, and to getting closer to the homolo- 
gous model. Furthermore, this decrease gets faster and 
faster, so that the central potential becomes — 00 after a 
finite timelapse: the cluster reaches the structure of the 
homologous model at a given time T = T2, and keeps it 
afterwards. 




Figure 12: Evolution of the central potential. 



Substituting (7.15) into (7.9), we obtain the explicit 
expressions of the fluxes near the center: 



dM 
dT 
<9H 
dT 



0.6464 
"U + 6.329 
0.6464 (U + : 
U + 6.329 



(7.17) 



The flux of mass is first positive and then decreases; it 
becomes zero at the time T 2 and remains null afterwards. 
The flux of energy is always negative; its absolute value is 
first decreasing; at T2, it reaches the value of -0.6464, that 
its keeps afterwards (see 5.22a). These results confirm the 
reasoning made in Chapter IV. 

At the beginning of its life, the cluster is likely not 
very concentrated; let's assume that Uo has initially the 
largest value allowed by (7.15), i.e. Uo = —6.329. By 
assuming that T = at the initial time, we find 



T2 is, in canonical variables 
central singularity. 



0.4427. (7.18) 
the time needed to form the 
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7.2 Proper differences 

We now suppose that the central density has already be- 
come infinite, but that the cluster as a whole has not 
completely reached the final state represented by the ho- 
mologous model. 

Let F be the distribution function of the homologous 
model; the distribution function of the cluster considered 
here would be 

F = F + AF, (7.19) 

where AF is small with respect to F . AF is time- 
dependent. In the same way, we set 



Q = Qo + AQ, etc. 



(7.20) 



The function F is assumed to be in the usual canonical 
form, defined by (4.6). The parameters K, b, c, would all 
have slightly different values than those of the homologous 
model: 

K = K + AK, etc. (7.21) 

The cluster obeys the relations (5.2) and (5.3), as soon 
as the terms including time derivatives, which were zero 
for the homologous model (see Equations 2.37e, 2.53, 3.12 
and 3.13), are put back in (5.2e) and (5.2f). These equa- 
tions become 



= 5" 



|»-|e)«J + (|c-f. 



(7.22) 



= (3A + 1)6+ (2 - 2A)c 

1A dln(M e ) d din (Re 



dT 



dT 



By substituting the expressions (7.19) to (7.21) into 
the previous equations, and neglecting the second order 
terms, we obtain a set of linear equations that have the 
dimensions of AF, AQ, etc., and their derivatives with re- 
spect to time. We know that the general solution of such 
a system is (unless a degeneracy exists) a linear combina- 
tion of particular solutions of the form 



-sT 



AF = AF e 
AQ = AQ c~ sT , 



(7.23) 



etc, 



where AF , AQ , ... do not depend on the time, and s is 
a constant, real or complex. As in the classical terminol- 
ogy, we shall call such a solution proper difference, and s, 
the associated proper value [T.N.: eigenvalue]. By substi- 
tuting (7.23) into the equations, and dividing by e~ sT , we 
obtain a system that is independent of time, but which 
involves a new parameter: s. 

7.2.1 Method of solution 

In practice, instead of writing and solving the system of 
equations of AF and AQ, etc it is easier to keep the 
equations of F, Q, etc to compute a slightly different 



model from the homologous one by means of these equa- 
tions, and to obtain the differences by simply computing 
the differences between the values of the two models. The 
differences should not be too large (so that the second or- 
der terms arc indeed negligible), nor too small (so that 
one can get them with a sufficient accuracy); experiment 
led us to fix the amplitude of these differences by setting 



AK = -0.001. 



(7.24) 



Taking into account that the device we use computes 
with 8 significant digits, we can obtain the differences 
with a precision of the order of 1/1000. 

Hence, the system to solve is as (5.2) and (5.3), as 
soon as (5.2e) and (5.2f) are replaced with 



(7.25) 




-(^-^]FQ 



FiQi dEi 



+s J [Qi(Fi - F ) - Fi(Qi - Qo)] dEij 



= (3A + 1)6+ (2 - 2A)c 



(7.26) 



+4s 



Re — Re 

Rp.n 



- A- 



IVL-IVL 



and when (7.24) is added. 

The method of solution is an extension of the one used 
for the homologous model (Chapter V): 

1. choose a value for s; 

2. choose a temporary form for the function F that sat- 
isfies (5.3a) and (5.3c), with K set by (7.21) and 

(7.24) ; 

3. compute D, Z, R, Q; 

4. adopt temporary value for 6 and c, and integrate 

(7.25) ; 

5. re-do after changing 6 and c, until the final conditions 
(5.3a) and (5.3b) are fulfilled; 

6. go back to point 3 with the new function F; 

7. when two consecutive approximations of F are equal: 
compute the right-hand side term of (7.26). It is not 
zero in general. Modify the value of s, go back to 
point 2; grope around this way with s until (7.26) is 
true. 

The computation is quite long, because of the trial and 
error required on three parameters: 6, c, s; that is why we 
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limited ourselves to the exact computation of the proper 
differences corresponding to the two smallest proper val- 
ues s. (We will see later that the possible values of s are 
real, positive and form a discrete series.) These differ- 
ences are the most interesting in practice, because they 
are those which decay the slowest. The other differences 
will be treated in a more approximate way later in this 
Chapter. 

7.2.2 Results: first proper difference 
We find that the smallest proper value is: 

s = 1.81. (7.27) 

The differences AF, AD, AR, AQ, AD P of the four 
fundamental functions and the projected density are given 
in Table |4] The differences of the parameters and the 
characteristic quantities are: 



Ab 



AR e 

AK 



= -6.62 



AK 
AM 



= +1.30 



AK 



= +5.22, 



(7.28) 



AK 



= -0.537 



AL e 

AK 



= +0.065. 



Figure Q~3] compares the projected densities of the 
usual homologous model (AK = 0) with those of the ho- 
mologous model modified by the first proper difference; 
here we have set AK = 0.2 so that the difference is well- 
visible. (As for its amplitude, the sign of AK is arbitrary; 
we could have chosen a negative AK; in this case, the dif- 
ference with the homologous model would have been in 
the opposite sense.) We see that the difference mostly 
affects the external regions of the cluster: the first proper 
difference mainly consists in a variation of the external 
radius of the cluster. The structure of the internal region 
differs little; the total mass and kinetic energy vary much 
less than the radius, as shown by the values (7.28). 

This is naturally explained by the fact that the pertur- 
bations between the stars are more efficient in the center 
of the cluster, where the density is higher; as a conse- 
quence this region becomes close to the final state, repre- 
sented by the homologous model, earlier. 

7.2.3 Results: second proper difference 
The second proper value is: 

s = 5.12. (7.29) 



We only give the differences of the parameters: 



Ab 
~AK 



= +1.71 



AK 



= +6.90, 



(7.30) 



AR P AM e AL P 

= -3.48 -r^- = -0.760 -r-^r = +0.182. 



AK 



AK 



AK 



The profile of the projected density (not shown here) 
shows that the second proper difference, like the first one, 
mostly consists in a variation of the external radius. 




Figure 13: Usual homologous model (AK = 0) and homologous 
model modified by the first proper difference (AK = 0.2). 



7.3 Stability of the homologous model 

Up to now, we have implicitly assumed that the clusters 
naturally tend toward the homologous model. For this to 
be true, the model has to be stable, i.e. any difference 
with respect to this model is decreasing. On the other 
hand, the stability of the homologous model, if demon- 
strated, will be, if not a rigorous proof, at least a very 
strong clue that the model is indeed the final state to- 
ward which all clusters tend. 

For the homologous model to be stable, it is necessary 
and sufficient that all the proper values s yield a positive 
real part. Thus, we are going to study the complete fam- 
ily of the proper values s, thanks to a very approximate 
calculation, yet sufficient for the goal we seek. 

We first assume that the central expansion (4.34) of F 
is valid up to E = —1. (This approximation, as all which 
will follow, has been suggested and checked by the exact 
computation of the first two proper functions.) Hence, 
for E = — 1, we have: 



AF , N 
AK^ 



AF' 
~AK 



(-1) 



J/2 



(7.31) 



Between E = — 1 and E = 0, we consider the ex- 
act equation (7.25). By calculating the derivative of this 
equation, neglecting the terms in F^" 1 ) and F 2 (because 
F vanishes at the boundary), replacing (7.19) etc and 
subtracting the equation of the homologous model, we 
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Table 4: 



"171 T T 

Ej Oi U 


AF 


AD 


AR 


AQ 


AD P 


AK 


AK 


AK 


AK 


AK 


-5 


12.18 


33.4 


0.00673 


0.0000287 


8.84 


-4.9 


11.62 


31.2 


748 


361 




-4.8 


11.08 


29.4 


831 


448 


8.42 


-4.7 


10.55 


27.6 


923 


565 




-4.6 


10.04 


25.8 


0.0102 


703 


7.98 


-4.5 


9.56 


24.2 


114 


885 




-4.4 


9.10 


22.6 


126 


0.000110 


7.54 


-4.3 


8.66 


21.1 


140 


138 




-4.2 


8.24 


19.7 


155 


171 


7.10 


-4.1 


7.84 


18.4 


172 


213 




-4 


7.45 


17.2 


191 


265 


6.66 


-3.9 


7.08 


16.0 


212 


329 




-3.8 


6.73 


14.8 


235 


408 


6.21 


-3.7 


6.40 


13.8 


261 


508 




-3.6 


6.07 


12.8 


290 


629 


5.77 


-3.5 


5.77 


11.9 


323 


782 




-3.4 


5.47 


11.0 


359 


970 


5.33 


-3.3 


5.20 


10.1 


399 


0.00121 




-3.2 


4.93 


9.33 


445 


150 


4.89 


-3.1 


4.68 


8.59 


495 


186 




-3 


4.44 


7.89 


552 


232 


4.45 


-2.9 


4.21 


7.23 


616 


289 




-2.8 


3.99 


6.61 


688 


360 


4.02 


-2.7 


3.78 


6.03 


769 


449 




-2.6 


3.58 


5.49 


861 


562 


3.59 


-2.5 


3.39 


4.98 


965 


703 




-2.4 


3.21 


4.50 


0.108 


882 


3.18 


-2.3 


3.04 


4.06 


122 


0.0111 




-2.2 


2.87 


3.64 


137 


139 


2.77 


-2.1 


2.71 


3.25 


154 


176 




-2 


2.56 


2.90 


174 


223 


2.37 


-1.9 


2.42 


2.56 


197 


282 




-1.8 


2.28 


2.25 


223 


359 


1.98 


-1.7 


2.15 


1.97 


253 


458 




-1.6 


2.02 


1.71 


288 


587 


1.62 


-1.5 


1.90 


1.47 


329 


755 




-1.4 


1.78 


1.25 


377 


975 


1.27 


-1.3 


1.66 


1.06 


433 


0.127 




-1.2 


1.55 


0.877 


500 


165 


0.95 


-1.1 


1.44 


0.718 


580 


217 




-1 


1.33 


0.576 


675 


287 


0.665 


-0.9 


1.22 


0.452 


791 


382 




-0.8 


1.11 


0.345 


933 


514 


0.418 


-0.7 


1.00 


0.253 


1.11 


700 




-0.6 


0.886 


0.177 


1.34 


965 


0.220 


-0.5 


0.769 


0.116 


1.63 


1.35 




-0.4 


0.646 


0.0684 


2.01 


1.93 


0.0816 


-0.3 


0.512 


0.0346 


2.55 


2.83 




-0.2 


0.364 


0.0130 


3.32 


4.28 


0.0111 


-0.1 


0.195 


0.0024 


4.53 


6.80 













6.62 


11.61 
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obtain 



FiQidEi • AF'' 



(7.32) 



+F"A ( f F1Q1 dEi J + I FiQidEi-AF' 

\J — oo I J— oo 

+F ' a (/-' 



F1Q1 dEi 



+ F'Qi(|l»-|c) + (| ( ,-|c)QAF' 
+ f h - 3 -c) F'AQ - FQ'Ac - cQ'AF 



cFAQ' + sQ'AF - sF'AQ j 

[T.N.: The upper limit of the first integral is missing in 
the original version] 

Furthermore (7.26) becomes, when neglecting the 
term in AM e : 



= (3A + 1)A6 + (2 - 2A) Ac + 4s 



AR e 



(7.33) 



The expansions (4.35) show that we have, near the 
center, 



R V2 



(7.34) 



AQ 

AQ' 



9^3 
3V3 



e E / 2 AK, 
e E / 2 AK. 



We shall suppose that these expressions are valid up 
to the boundary. In (7.32), we can neglect the second 
and the fourth terms, because the integrals do not change 
much. We also neglect the terms in Ab, as they are small 
with respect to Ac. By taking Ac from (7.33) and sub- 
stituting it in (7.32), we obtain 



| J FiQi dEi • AF" 



(7.35) 



FiQi dEi 



+(s - c)QAF + (s - c) 
-^(fQ^F'Q 



3A + 1 
9^ 



c Q AF' 



AK 



e E / 2 F'Q 



+c 



6^+4 9^3^/2 F'Q + ^2e E / 2 FQ' 
3A + 18V2 2^2 



t^(fq' + ^ q 



AK 



Finally, we neglect the term in AF', whose factor is 
relatively small, and the term in c. The equation shrinks 
to: 

= AF" + (s - c)B l AF + (s - c)B 2 AK, (7.36) 

where B\ and B 2 are two functions of E. This is a sec- 
ond order differential equation for AF; it must fulfill the 
boundary conditions (7.31), as well as the condition 



AF(0) = 0. 



(7.37) 



These three conditions would be simultaneously ful- 
filled only for certain values of s, which are the proper 
values we seek. 

The numerical computation of the functions B\ and 
B-i shows that they can be quite well fitted, in the range 
(-1,0), with the expressions 



By 

B 2 



2.10 



(0.332- E) 2 ' 
(0.265 + 0.643 E) B x . 



(7.38) 



By replacing this in (7.36), we find that the general 
solution of this equation is 



AF 

~AK 



-0.265 - 0.643 E (7.39) 
+Ci(0.332 -E) Pl + C 2 (0.332 -Ef 2 , 



where C\ and C 2 are two arbitrary constants, and pi, p 2 
are the roots of the equation: 



p 2 -p + 2.W(s - c) = 0, 



so that we have 



P1+P2 = 1, 
Pl P2 = 2.10(s - c) 



(7.40) 



(7.41) 



By writing the boundary conditions (7.31) and (7.37) 
for the formula (7.39), and then eliminating C\ and C 2l 
we obtain 



(1.332)p i (1.332)P 2 1.271 
Pi(1.332)p i p 2 (l-332)P 2 -0.181 
(0.332)^ (0.332)P 2 0.265 



= 0. 



(7.42) 



The equations (7.41a) and (7.42) make a system of 
equations for p\ and p 2 that can be solved numerically. 
We do not give the details of this solution, which is long 
but without difficulty; we find that p\ and p 2 are neces- 
sarily of the form 



Pi = 2 



(7.43) 



where £ is real, and i = y/—T. The possible values of 
£ constitute an infinite series; the first ones are given in 
the table below, as well as the corresponding values of s, 
found from (7.41b). 
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s 


1 91 


1 99 
l.ZZ 


o.oZ 


7 47 


5.43 


14.6 


8.28 


33.2 


9.91 


47.3 


12.8 


78.3 


14.4 


99.6 


17.3 


143 


18.9 


171 


21.8 


227 



That is, all the proper values s are real and positive; 
therefore, the homologous model is stable, and it is likely 
the final state toward which all the clusters tend. 

The fact that s never yields an imaginary part shows 
that the clusters tend toward the homologous model 
through a simple "relaxation" , with no oscillations. 

The first two values of s are in rough agreement with 
the exact values (7.27) and (7.29); the differences are not 
surprising, given the great number of approximations we 
made. 

Note that the values of s increase very fast which 
shows indeed that only the first values are interesting, 
in practice; the next ones corresponds to differences that 
vanish very quickly. 

Finally, we note that these results allow us to define 
rigorously the "relaxation time" for the entire cluster: this 
would be the time required for the amplitude of the first 
proper difference to be divided by e. 

7.4 General picture of the evolution 

We have identified three distinct evolutionary phenom- 
ena: 

1. the normal homologous evolution; 

2. the formation of the central singularity; 

3. the decrease of the differences. 

We are now going to compare their timescales. 

The normal homologous evolution can be illustrated 
through the variation of the total mass given by (5.15) 
and plotted as a straight line on Figure [TU The cluster 
disappears after a time T\ given in (5.13): 



Ti = 



1 

7o c 



2.452 

7o 



(7.44) 



The formation of the central singularity can be seen 
as the variation of AM, difference in mass with respect to 
the homologous model. This variation is given in (7.7a) 
and (7.16), in homologous variables; we switch to non- 
homologous variables by means of (5.12) and (5.15). The 
resulting curve is plotted in Figure[T4](for all the curves of 
this figure, the mass or the mass difference is normalized 
with respect to its own value at T = 0). 




Figure 14: Compared evolutions: EH = homologous evolution 
[T.N.: evolution homologique]; EP1 = first proper difference; EP2 
= second proper difference [T.N.: ecart propre]; C = formation of 
the central singularity. 



We immediately see that the formation of the cen- 
tral singularity is relatively rapid: it ends at the time T2, 
computed from (7.18) and (5.12): 



To 



0.4051 

7o 



(7.45) 



and thus occurs for only the first sixth of the lifetime of 
the cluster. 

Finally, the decrease of the proper differences will also 
be represented as variations of the mass differences. In 
homologous variables, we have, from (7.23), 



AM e 

AM,., 



-sT 



(7.46) 



and we switch back to non-homologous variables by 
means of (5.12) and (5.15), which gives 



AM e 
AM, 



1 



T 



l + s/c 



(7.47) 



This variation is plotted in Figure [T4j for the first two 
proper differences. We see once more, that the decrease 
of the differences is rapid with respect to the homologous 
evolution. 

In reality, the three evolutionary phenomena overlap. 
The initial amplitude of the differences with respect to the 
homologous model cannot be determined in the frame- 
work of this theory; it could only be found thanks to a 
detailed study of the formation and the initial phases of 
the evolution (see Chapter IX). However, there is no rea- 
son for the initial state of the cluster to already be close 
to the homologous model, and thus the initial differences 
in the various quantities are likely of the same order of 
magnitude as the quantities themselves. Therefore, we 
can suppose that the curves of Figure 1141 normalized to 
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unity for T = 0, roughly represent the relative importance 
of the various effects. 

The evolution of a cluster can be describe in a broad 
outline as follows: the central density first increases, 
faster and faster; simultaneously, the differences with re- 
spect to the homologous model decrease everywhere in the 
cluster. After a time Ji/6, the central density has became 
infinite; the remaining difference is practically limited to 
the first proper difference; the cluster has only lost 1/6 
of its mass, through escapes. Then, the first proper dif- 
ference continues to decrease; after a time Ti/3, it has 
almost vanished; the cluster takes the form of the homol- 
ogous model and keeps it until it disappears, at the time 
T x . 

We will see in the next Chapter how these results can 
be used to estimate the age of globular clusters. 

8 APPLICATION TO GLOBULAR CLUS- 
TERS 




Figure 17: Spacial density of the cluster M3, from Kholopov (top) 
and Oort & van Herk (bottom). The curves are vertically shifted. 



We are going to compare the theoretical results from the 
previous Chapters to observational data of globular clus- 
ters. The first goal of this comparison is to check that the 
theory is compatible with the facts (but this would not 
be a proof of the correctness of the theory, for the rea- 
sons detailed in Chapter I); the comparison will give us 
information that is out of reach of the direct observation, 
in particular about the mass-luminosity relation and the 
age of the clusters. 

8.1 Artificial cluster 

In order to get a concrete picture of the theoretical model, 
an "artificial cluster" has been built by calculating the 
spacial coordinates of the stars, from a table of random 
numbers, so that it reproduces the density law of the ho- 
mologous model. More precisely, let m, n 2 , n3 be three 
random number, taken between and 1, with an uniform 
probability density function; the spherical coordinates R, 
6, <p of the star are computed by means of 



The central conde nsation in Figure [TH] is likely to be too 
high (|Kindll961ah . 



M R = mM e , 
cos 9 — 2n-2 — 1, 
ip = 2irri3 . 



(8.1) 



Figure [15] plots the two-dimensional projection of the 
cluster obtained. The number of stars is 1320. The figure 
covers only a fraction of the surface of the cluster (see 
the scale at the bottom of the figure), but, however, con- 
tains all the stars; this is because the projected density is 
extremely small in the external regions. 

Figure [TC] is a zoom-in on the central region of the 
cluster. 

Globally, this artificial cluster resembles well the real 
clusters; the similarity would not be perfect because the 
artificial cluster of Figure [15] corresponds to the case of 
equal masses, while the real clusters yields an extended 
mass spectrum, of which we only observed the upper end. 



8.2 Comparisons of the projected densities 

A more precise comparison of the homologous model with 
real clusters can be obtained when plotting the density 
profiles. It seems better to do the comparison for the 
projected densities, and not the spacial densities. Indeed, 
the theoretical curves can be obtained with an arbitrary 
accuracy and going from the theoretical spacial density 
to the theoretical projected density is errorless; on the 
contrary, going from the observed projected density to 
the spacial density strongly amplifies the errors, that were 
already not negligible. This is shown in Figure [T71 where 
the spacial density of the cluster M3, derived from the 



observations thanks to two d ifferent methods (jKholopov 
19551: lOort fc van Herk]|l959h is plotted. 

Generally, it seems better to transform the observa- 
tional data as little as possible, therefore to compare the- 
ory and observations not at an intermediate level, but 
rather at the very level of the observations. 

We shall use, not the projected density itself, but its 
product with the distance to the center: ppr. This has 
several advantages: 

1. in the theoretical model (case of equal mass stars), 
this quantity tends to a finite value in the center, 
while the projected density becomes infinite. 

2. ppr varies more slowly than the projected density 
(recall Figure [5] and Figure [9]); 

3. the observations are often counts of stars in concen- 
tric rings of constant width. The number of stars in 
a ring is 

1 



An = 



2npp r dr. 



2) 
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Figure 15: Homologous model: artificial cluster. 
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Figure 16: Zoom-in on the center of Figure IT5l 
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Because this function is smooth and does not yield a 
strong curvature, we can write 

An = — Ar (p P r)-, (8.3) 



rr> 



with: 



Ar = r 2 -n, r = (8.4) 



That is, the observed numbers An immediately give, 
up to a factor, the values of ppr. 



Table 5: 



M v 


My 


My 


/' 


log (CO 


MA) 


M bo] 


min 


max 










— CXI 


0.4 


-0.5 


1.55 


1.328 


0.090 


-0.6 


— CXI 


1.2 


-0.1 


1.50 


1.538 


0.076 


-0.2 


— cx> 


1.9 


0.5 


1.45 


1.698 


0.061 


0.4 


— cx 


3.5 


2.2 


1.42 


2.472 


0.052 


2.1 


— cx 


4.1 


2.8 


1.40 


2.672 


0.046 


2.7 


— cx 


1.6 


3.3 


1.10 


2.894 


1.941 


3.2 


— cx 


5.4 


.3.9 


1.05 


3.125 


1.921 


3.8 


— cx 


6.3 


4.5 


0.80 


3.246 


1.803 


1.4 


4.6 


6.3 


5.5 


0.60 


2.986 


1.678 


5.4 



8.2.1 Projected density: star counts (MS) 

Sandagd (|l954lll957al) has made detaded counts of stars 
in the cluste r M3 = NGC 52 7 2; th e results have been 
published bv lOort fc van Herldll 1959ft . These counts cor- 
responds to the stars brighter than a given magnitude, 
and situated in concentric rings. It would be necessary, 
in principle, to derive the numbers of stars in succes- 
sive magnitude ranges by subtraction; however, attempts 
showed that this strongly increases the spreading of the 
points, because the errors accumulate, while the number 
decreases. Furthermore, because the lu minosity functio n 
rapidly grows with the magnitude (see Sandaee 1957aft . 
the stars brighter than a given magnitude are almost all 
gathered in a small range of magnitude, and thus all have 
similar masses (we assume a well-defined mass-luminosity 
relation) . It is therefore better and allowed to accept that 
the observed numbers correspond, for each limit in mag- 
nitude, to stars of a given mass. 

That is, the points on Figure [18] mark the observed 
numbers of stars brighter than a given magnitude My, 
in rings of width 30" situated at the mean distance r 
from the center. These numbers must be compared to 
the theoretical curves obtained in Chapter VI (Figure . 
By trial and error, we find that the best match is obtained 
for 

V ~ 0.21 canonical units. (8-5) 
We note that this gives out an external radius 

r e = 4.703 canonical units ~ 22.4' ~ 79 pc. (8.6) 

For each group, the curve is fitted, by playing with two 
parameters: the relative mass and factor C, i.e. the 
number of stars in the group; a variation of C corresponds 
simply to a vertical shift of the curve. 

The values of p and C determined this way are given 
in Table columns 4 and 5. The curves are plotted in 
Figure [TU we see that the agreement is as good as pos- 
sible, given the accuracy of the observations. Only the 
brightest stars are not well-retrieved near the center of 
the cluster, which can easily be explained by two reasons: 
(1) we have seen in Chapter VI that the simplified theo- 
retical model is likely to be wrong near the center, where 
it predicts a too high density of massive stars; (2) the ob- 
servation underestimates the number of stars in the very 
populated central regions. 




Figure 20: Mass-luminosity relation in M3, from the spacial dis- 
tribution of the stars (dots) and from the theory of stellar evolution 
(line). 



8.2.2 Mass-luminosity relation of the stars 

Table [5J column 3, gives the mean absolute magnitude 
My of each group, c omputed from the luminosity func- 
tion (|Sandagdll957al ). We notice that the difference be- 
tween this mean magnitude and the maximum magnitude 
is never big, which confirms what has been said above. 

The values of p and My allow us to plot the mass- 
luminosity relation of the cluster (Figure [2"0"|) ; p is propor- 
tional to the mass m of the stars, according to (6.32). The 
bolometric correction has been uniformly taken equal to - 
0.1. In order to extend the relation as much as possible for 
the faintest stars, we made an exception and considered 
the group of stars whose magnitude is between 4.6 and 6.3 
(in the case of the narrower group 5.4 < My < 6.3, the 
dispersion of the points is too large and we cannot plot a 
curve). This group is compared to the closest theoretical 
curve in Figure 1191 the parameters are given in the last 
row of Tabled 

The points in Figure |2"U1 show a well-defined relatiorH. 
We observe a strong bend at My = 3, which corresponds 
exactly to the bend in the Hertzsprung-Russell diagram, 
and which marks the beginning of the zone of rapid evo- 
lution of the stars; the stars with a magnitude below 3 



^von Hoernerl 1^957, Figure 8) has obtained a similar curve, from 
the same observations, but with a different theoretical model. Our 
results indicate a larger variation of the mass as a function of the 
luminosity. 
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Figure 19: Continuation of Figure ITSl 



almost have all the same mass. The slight increase of the 
points on the left-hand side is likely not real and comes 
from errors when computing [i. 

It is interesting to compare this mass-luminosity rela- 
tion, deduced from the purely mechanical consid erations , 
to those provided by the stellar evolution theory. Sandaeel 
(|l957bl Table 6) has computed the relation between the 
present-day magnitude and the initial magnitude of the 
stars of M3, thanks to a semi-empirical method based on 
the observed HR diagram and on the theoretical evolu- 
tionary tracks. We get the mass from the initial magni- 
tude, assuming that the initial mass-luminosity relation 
is that of the main sequence, i.e. 



find that 



L 



Mr. 



•7) 



This way, we get the present-day mass-luminosity re- 
lation, the solid line on Figure [501 the two relations 
have been vertically shifted so that their horizontal parts 
match. 

The two curves are clearly in disagreement. The posi- 
tion of the bend is more or less the same; but on the low 
luminosities end, the "dynamical" masses decrease much 
faster than the mass found from the evolution. 

It is hard to tell which curve is wrong. Neither is very 
reliable. The finding of the dynamical masses relies on an 
approximate calculation (Chapter VI), which may lead to 
large errors when (as is the case for globular clusters) the 
mass spectrum is very broad. On the other hand, the hy- 
pothesis stating that the initial mass-luminosity relation 
of the cluster would match the main sequence, should be 
taken carefully, because of different chemical abundances. 

We will explore the possibilities of progress in the com- 
putation of the dynamical masses in the last Chapter. 

The adjustment of the vertical scales gives an inter- 
esting piece of information: when compared to (6.32), we 



— = 0.80 M 



©• 



This value can also be computed from the mass spec- 
trum of the clust er. From the ini tial luminosity function 
ijj(M v ) given bv lSandagel (|l957al Ta ble 2 an d Figu re 2) 
and the mass-luminosity relation of iKuiperl (Il942h . we 
obtain the initial mass spectrum. The present-day spec- 
trum is derived by assuming that all the stars heavier than 
1.44 Mq have been evolving down to this value through 
their transformation into white dwarves. The effect of 
escape can be neglected. We find 



m = 0.353 Mq, 
7^2 = 0.265 M , 
Tr^/m = 0.75 M Q , 



(8.9) 



in very good agreement with the value (8.8). 

This previous finding is quite sensitive to the mass 
chosen for the white dwarves, TOf, [T.N.: The subscript 6 
stands for naine blanche: white dwarf.]. For m^ = 1 Mq, 
we find to 2 /to = 0.61 M Q ; for to;, = 2 Mq, we find 
to 2 /to = 0.98 Mq. The comparison of the two findings 
seems to confirm that the white dwarves have the limit 
mass of Chandrasekhar, mi, = 1.44 M Q . However, here 
again, the uncertainties of the theory and of the determi- 
nation of the mass spectrum question the validity of this 
conclusion. 

8.2.3 Projected density: brightness measurements 
(47 Tuc) 

iGascoigne fc Burn ( 1956 ) have obtained the curves of pro- 
jected density of the clusters 47 Tucanae and u> Centauri 
from a photometric method. The disadvantage of this 
procedure is to provide only the total brightness, without 
distinguishing the stars of different magnitudes; however, 
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Figure 21: 47 Tuc: comparison between observations (dots) and 
theory (line). 



it allows for a precise measurement of the density from 
the center of the cluster up to quite remote distance. 

The case of u Centauri will be studied later. Figurel2"Tl 
plots the observed values for 47 Tuc = NGC 104 and the 
closest theoretical curve, which corresponds to: 

1' = 0.088 canonical units, (8.10) 
/i = 1.33, 

C = 0.555 stars of the tenth magnitude per minute. 
The external radius is then: 



53.4'. 



(8.11) 



The value of /i found here can been seen as a partic- 
ular "mean mass" , obtained by weighting the stars pro- 
portionally to their luminosity. It is only slightly lighter 
than the mass of the heaviest stars (fj, ~ 1.45 from the 
results obtained for M3) which contribute to the major 
part of the total brightness. 

Figure [5T] shows good agreement (the small systematic 
deviations of the points from the theoretical curve can be 
explained: the data of Gascoigne & Burr does not directly 
come from observations, but has rather been read on a 
smoothed curve). Near the center, the theoretical curve 
is too high, for reasons already listed in relation to M3. 

Gascoigne & Burr have measured the total luminosity 
in a series of circles of small radius around the center of 
the cluster. For each circle, we can compute the mean 
value of ppr, by using (8.3), with r\ = 0; the crosses on 
Figure [5T] mark the value obtained. We clearly see that 



ppr tends toward a well-defined central value. This seems 
to be a good confirmation of the shape predicted by the 
theory near the center, and in particular of the fact that 
the central density is infinity 

8.3 Mass-radius relation 

The relation (3.7) links the mass of a cluster with its ra- 
dius and the mean value of the negative curvature of the 
galactic field along its orbit. Unfortunately, the exact mo- 
tions of the globular clusters in the Galaxy are unknown; 
the r adial velocities can be mea sured with a good preci- 
sion (|MavalJll94ilKinmanlll959l ). but the proper motions 
stand at the limit of the observational possibilities, an d 
have been obtained for 9 clusters only (jGamaleil Il948l ). 
with a low precision. It seems even not possible to know 
wether the orbit s of the clusters ar e in average rather ra- 
dial or circular ( von Hoerner 1955 : Kurth 1960f) . That is 
why we will consider the curvature of the galactic field 
to be the same for all the clusters, for want of anything 
better. 

The radius and the mass of the globular clusters 
should be linked, at least approximately, through a re- 
lation like (3.8). We are going to check this using obser- 
vational data. To limit the effect of observational errors, 
we will only use a series of homogeneous observations, 
done by the same author, under the sam e circum s tances . 

The distances have been taken from iKinmanl (|1958l ). 
In order of preference, we take the distances derived from 
the position of the main sequence, from the magnitude 
of the variable stars, or from the magnitude of the 25 
brightest stars. If none of these three estimates exists, 
the cluster is not considered. 

Christie ( 1940l) made very precise measurements of 



the apparent magnitudes of the clusters. From them, we 
derive the absolute photographic magnitudes M pg ; the 
possible effect of absorption is cancelled, because the dis- 
tances themselves are derived from the apparent magni- 
tudes of the stars. 

Then, we assume that the ratio of the total mass to the 
total luminosity is the same for all t he globular clust ers; 
by taking the mass estimate of M3 bv iSandagd (|1957al ) as 
reference, we find that the mass A4 e of a cluster is given 
as a function of its absolute magnitude M pg by 



log .M e 



-0.4 M v 



2.00. 



.12) 



This is confirmed by the study of Zeliakhl (1957) who, 
after a detailed discussion on several clusters, arrived at a 
quite similar value for the constant: 1.94 instead of 2.00. 
Finally, we take the app arent external radii measured 



r many, we take t nc apparent external rami measured 
by IShaplev fc Saverl ( 19351 ): they appeared to be more 
precis e that those published more recently by [Mowbray 
(|1946l ). To derive the real radii, we should, in principle, 
use the distances corrected for absorption. However 



as 



9 However, rKinel jl961bT) has observed the central region of sev- 
eral clusters in detail and concluded that a high but finite central 
density exists. 



38 



M. Henon, translated by F. Renaud 



NGC 



Table 6: 

log (Me) 



log (re) 



288 


4.64 


1.40 


1904 


5.12 


1.37 


2419 


5 08 


1.71 


4147 


4.44 


1.03 


5024 


5.12 


1.66 


5139 


5.80 


1.90 


5272 


5.40 


1.64 


5897 


4.60 


1.50 


5904 


5.28 


1.60 


6093 


5.04 


1.52 


6121 


4.48 


1.30 


6171 


4.36 


1.34 


6205 


5.20 


1.38 


6218 


4.80 


1.48 


6229 


4.88 


1.38 


6254 


5.04 


1.53 


6266 


5 28 


1.60 


6273 


5 04 


1.50 


6284 


4.56 


1.49 


6293 


4.84 


1.39 


6333 


5 00 


1.53 


6356 


5.16 


1.70 




o.uu 


j-.oy 


6626 


4.84 


1.46 


6638 


4.88 


1.32 


6656 


5.08 


1.54 


6715 


5.36 


1.71 


6723 


5.00 


1.30 


6779 


4.68 


1.38 


6809 


4.84 


1.46 


6864 


5.48 


1.74 


6981 


4.64 


1.42 


7078 


5.44 


1.60 


7089 


5.56 


1.63 


7099 


4.76 


1.29 



indicated bv lShaplev k Saverl (|l935h ; IShaplevl (|l949h . ab- 
sorption induces two opposite effects: an overestimate of 
the distance and an underestimate of the apparent radius. 
A quantitative study show s that these two effects balance 



each other almost exa ctly ([Parenago. Kukarkin. k Floria 



19491 lLohmmm1ll952h . Therefore, we can completely ne- 
glect the possible existence of absorption. For that mat- 
ter, this is very fortunate because the absorption is not 
quite well known. 

Table [5] gives the list of the 35 clusters for which 
the three measurements exist, with their mass (in solar 
masses), and their radius (in parsecs). The mass-radius 
diagram is plotted in Figure [^U we see that, despite the 
dispersion of the points, a defined relation arises. The 
main source of error comes from the estimation of the 
external radius, which is difficult to observe precisely; 
this error alone is enough to explain the dispersion of 
the points, and it is even remarkable that it does not in- 
duce a larger dispersion. (This shows the quality of the 
meas urements of Sh apley k Sayer; if one uses the radii 
from lMowbravlll946l the dispersion is much larger.) 

In addition, the small dispersion justifies the hypoth- 
esis of constant curvature of the galactic field for all the 
clusters. 

The largest error being on r e , we shall compute the 




1,6 1,8 

logr e 



Figure 22: Mass-radius relation of the globular clusters. 



means in the horizontal direction. This way, we find that 
the data points are well fitted with the relation 



log (r e ) = ~ log (M e )- 0.17 (±0.11) 



i.e. 



3.2 M Q /pc 3 



(8.13) 



(8.14) 



in agreement with the theoretical formula (3.8). 

The isolated point in the bottom-left corner of the fig- 
ure is NGC 4147. This cluster is easily distinguishable 
from the others; Shapley k Sayer qualify it as "abnor- 
mal" . The point on the top-right corner is uj Cen, also 
abnormal, as we will see later. 

8.3.1 Correction of the external radii 

The observed external radius is, in fact, systematically 
smaller than the real radius. Indeed, the projected den- 
sity decreases very rapidly toward the outskirts and be- 
comes very small long before reaching the boundary of 
the cluster. This clearly appears in the artificial cluster 
of Figure [T5J a direct estimate of the radius of this clus- 
ter, from its aspect, gives a value of the order of 0.6 times 
the real radius, and this without any background effect. 
Table Q] indicates that only 3% of the total mass of the 
cluster is situated (in projection) beyond half of the ra- 
dius. That is, even without absorption, the radius may 
be underestimated by a ratio of about 1/2. A good con- 
firmation of this value can be obtained in the case of the 
two clusters studied above: for M3 and 47 Tuc, the ob- 
served radii are 11' and 28.2', i.e. 0.49 and 0.53 times the 
real radii derived from the theory and the projected den- 
sity profiles, and given in (8.6) and (8.11), respectively. 
Furthermore, there is no absorption for these clusters. 
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Thus, we will systematically make a correction, as- 
suming that the real radii are twice as large as the ob- 
served radii, given by Shapley & Sayer. As a consequence, 
the relation (8.14) must be replaced with 



— =2. = 0.40 M /pc 3 . 



(8.15) 



This relation between the mass and the real external ra- 
dius is plotted in Figure as a dashed line. All the 
observed radii stand on its left side, as expected. 

8.3.2 Curvature of the galactic field 

When comparing (3.7) and (8.15), we find the mean cur- 
vature of the galactic field for the clusters: 



(PUc 
dx 2 



-5.4 x 10~ 16 yr~ 2 



-520 kmV 2 kpc~ 



(8.16) 

For comparison, we can compute the curvature in the 
solar neighborhood; it derives from the Oort's constants A 
and B (the numerical values are taken from I Alien! 1 1 9 5 51 ) : 



d 2 u G 

dx 2 



(B - A)(3A + B) = -1370 kmV 2 kpc~ 



(8.17) 

The mean curvature of the galactic field for the glob- 
ular clusters is thus about 3 times smaller than in the 
solar neighborhood. This result is plausible, because the 
clusters travel in average at a quite long distance from 
the center of the Galaxy, and not in its plane. A more 
detailed comparison of the value (8.16) with the known 
structure of the Galaxy would be interesting, but is out 
of the scope of this work. 

8.4 Evolution 

Knowing the mass M e and the radius r e of a real cluster, 
we can compute the parameters j3 and 7 of the homol- 
ogy, using the transformation relations (6.31b), (6.34a), 
(2.32c), (2.43a). We get 



7 



G 



Mt 



(lWG -3, MS , 



M. 



-1/2 



(8.18) 

-3/2 



This allows us to go from any quantity of the homolo- 
gous model to the corresponding physical quantity, using 
the transformation formulae (2.31), (2.32), (2.43), (2.48), 
(6.31) and (6.34). In particular, we obtain: 



dt = r -i/2 m 1 
dT 



2 ln(n) \M 



a a \ 1/2 / \ 3/2 



R, 



(8.19) 



This relation gives the time scale of the evolution of 
the cluster, because T is, as we have seen (Chapter V), 
a parameter that measures the level of evolution of the 
cluster, while t is the real time. M e and R e are given in 



(5.5). We use the value to 2 /to found from the observation 
of M3 (Equation 8.8). The factor In (n) does not vary 
much for the different globular clusters; we use the typical 
value n = 10 5 . Finally, the radius is linked to the mass 
through (8.15), found from the observations. Thus, we 
find the very simple relation: 



H = 1.76 x 10 5 M c 



.20) 



where t is in years and M e in solar masses. 

(6.22) provides the escape rate of the stars of given 
mass, which also reads 



1 dn m 
^~dT 



.21) 



Recall that n m is the number of stars whose mass is be- 
tween m and m + dm. 9 is given in Table [5] and can be 
fitted with the relation: 



0- 0.4078 (3 - 2/i). 



.22) 



Then, we can compute the variation of the total mass, 
by using (6.25b), (6.32) and (6.26). We find 



dM e 
dT 



-0.4078 M t 



.23) 



We note that this result does not depend on the mass 
function. 

When comparing with (8.20), the mass Ai e cancels, 
and we get 



dM t 
dt 



-2.3 x 10" 6 M Q /yr. 



.24) 



Thus: the mass of a globular cluster decreases by 
2300 Mq per billion years; this rate is almost the same 
for all the globular clusters, and it also remains constant 
with time. The simplicity of this result is to be noted. 

In fact, the rate is not exactly constant. The term 
In (n) changes with time, which would lead to slight slow- 
down of the evolution near the end. On the other hand, 
the term too = m 2 /to increases because the lightest stars 
escape faster. This effect is opposite to the previous one 
and likely stronger. For example, in the cluster M3 which 
is still at the beginning of its evolution, the value observed 
is too = 0.80 M ; near the end of its evolution (in about 
100 Gyr), only the heaviest stars will remain and we will 
get mo = 1.44 M Q , thus a multiplication of the evolution 
rate by a factor 1.8. 

8.5 Mass function of the globular clusters 

The evolution law of the masses of the globular clusters 
partially sets the present-day mass distribution of these 
masses. We are going to see whether a confirmation of 
(8.24) can be obtained this way. 

From the apparent m agnitudes and th e apparent dis- 
tance moduli provided bv lLohmannl (1952) for 94 clusters, 
we compute the absolute magnitudes and then the masses 
using (8.12). 
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Figure 23: Mass function of the globular clusters, observed (dot- 
ted line), corrected for observational selection (solid line), initial 
(dashed line). [T.N.: The upper axis reads t en milliards d'annees: 
t in billion years.] 



The mass function we obtain is plotted in Figure [231 
using the dotted line on the left-hand side. The two most 
massive clusters are off the figure: M2 and uj Cen of 
masses 3.4 x 10 5 M and 5.2 x 10 5 M . We must also 
investigate the effect of observational selection. By look- 
ing at the apparent magnitudes m pg of the clusters, we 
see a clear upper limit at m pg — 12.5; we assume that 
this value marks the highest observable magnitude. As 
a consequence, the clusters of given absolute magnitude 
Mpg are observed only if their apparent distance modulus 
is smaller than 12.5 — M pg . That is, we can compute, for 
each value of M pg , the observable fraction of the clusters, 
as soon as we know the distribution of clusters as a func- 
tion of the apparent modulus. This distribution (that we 
assume to be independent of the absolute magnitude) is 
easily found thanks to the brightest clusters, for which 
observational selection does not play any role. 

The solid line in Figure [231 is the corrected mass func- 
tion. It can be fitted with a smooth curve. In particular, 
we note that this curve takes a finite value for j\4 e = 0. 
But if we assume that all the globular clusters formed 
in the past and that no more are born now; and if the 
evolution of the masses is as 



dMe 
d* 



.25) 



we can easily find that the present-day mass function 
must be proportional to M.~ v near the origin [T.N.: i.e. 
for j\4 e — > 0]. Therefore, this mass function must be zero 
for A4 e = if v is negative, and become infinite if v is 
positive; it can take a finite value only if v = 0. Therefore, 
the observed mass function seems to confirm the evolution 
law (8.24). 

The validity of this point is slightly weakened because 



it relies on the number of the least massive clusters, that 
are the least well observed; but we can argue that the 
correction for observational selection (Figure [23]) is not 
very important; the conclusion would have been the same 
without the correction. 

Note also that (8.25) is not fanciful but in fact covers 
all the evolutionary mo dels already studied by s everal au- 
thors. The formulae of IChandrasekhar (|l943bllc ) assume 
v = +1; the model from iKina (11958a) corresponds to 
v = —2.5 and the similar model of Ivon Hoernerl ([19581 ) to 
v = -3.55. 

The evolution formula (8.24) translates into a simple 
shift of the mass function toward the left; or alternatively, 
into a shift of the zero-mass toward the right. The upper 
axis of Figure [231 indicates the position of this zero-mass 
at different times, past and future. In particular, the 
[T.N.: vertical] dashed line on the left-hand side marks 
the position of the zero-mass point at the time of birth 
of the clusters, assuming it occurred 25 Gyr ago (see the 
next Section). By extrapolating the observed curve a lit- 
tle, we obtain the initial mass function. We find that 
about 62 globular clusters have disappeared since the ori- 
gin of the galaxy. Nowadays, one cluster disappears every 
500 Myr, in average. Every cluster that has survived has 
lost 57000 Mq as escaping stars since its birth. For all 
the clusters, the mass of the escapers since their birth is: 
8 x 10 6 Mq. 

In the future side, we find that half of today's clusters 
will be gone in 30 Gyr; but the most massive ones will 
survive much longer: uj Centauri will reach 230 Gyr. 

Up to now, we have negle cted the phenomeno n of gas 
ejection by the massive stars. Ivon Hoerner ( 1958 ) showed 
that this phenomenon is only relevant at the beginning of 
the life of the cluster: 30% of the initial mass is lost in 
the gaseous form within the first billion years, then only 
6% during the next 5 Gyr. The results obtained above 
are therefore valid, as soon as we put aside the initial 
evolutionary phase, lasting about one billion years. 

8.6 Age of the globular clusters 

When a cluster reaches the homologous state, any sign of 
its past (in particular its initial structure and its age) is 
lost. That is, one should specially care about the clusters 
that still show significant differences with the homologous 
model. According to the theory (see Equation 8.20), the 
most massive clusters are those which evolve the slowest, 
therefore, those where the largest differences should re- 
main. That is, we are interested in uj Centauri, the most 
massive of all globular clusters. 

But precisely, uj Centauri yields a peculiar structure, 
very different from those of the other clusters. The differ- 
ence already arises as its concentration class in Shapley's 
classification is VII, while all the other high luminosity 
clusters lie in classes I to V. The difference is even more 



visible in the work of iGascoigne fc Burn (jl956l ) who ob- 
served uj Cen and 47 Tuc under the same circumstances. 
The Table II and the Figure 2 from these authors per- 
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fectly show that the light from a small circle of radius r 
surrounding the center of the cluster is proportional to 
r 2 in the case of ui Cen, and to r for 47 Tuc. In other 
words: in 47 Tuc, the projected density goes as 1/r near 
the center; but in u Cen, it is almost constant. The cen- 
tral density of ui Cen is finite. 

The theoretical interpretation comes immediately: 
47 Tuc has already reached the stage in its evolution 
where the central density becomes infinite; uj Cen, be- 
ing more massive, has not reached it yet. (It seems to be 
the only globular cluster in this case.) 

This confirms the theory, and we now have a simple 
way of estimating the age of the clusters. When com- 
paring (2.32e), (7.44), (7.45) and (8.20), we find that the 
time required for the infinite central density to appear is 



U = 0.713 x 10 5 M t 



0.854 x 10 5 M e 



3.26) 



where ti is in years, M eo 1S the initial mass of the cluster 
and Me 2 hs mass at the time when the infinite central 
density appears. The present-day masses of 47 Tuc and 
uj Cen, obtained as in the previous Section, are: 2.8 x 
10 5 M Q and 5.2 x 10 5 M Q . We conclude that the age of 
these clusters (assuming it is the same) must be between 



24 x 10 9 



vi- 



and 



44 x 10 a 



yr. 



1.27) 



Furthermore, the most recent estimates of the age 
of the globular clusters from the stellar evolution theory 
(|Sandagelll96ll ) are 



22 x 10 9 yr for M13, (8.28) 
26 x 10 9 yr for M3 and M5. 

That is, two completely independent methods lead to 
values in agreement for the age of the globular clusters. 
This result is very encouraging, both for the theory of the 
dynamical evolution of the clusters and for the theory 
of stellar evoluti on. It also em phasizes the discrepancy 



underlined by ISandaed (|196lh between the ages of the 
globular clusters and the recent models of an expanding 
Universe. 

According to theory, u> Cen should have more differ- 
ences with respect to the homologous model; the first 
proper difference should still exist with a non negligible 
amplitude. We have seen (see Figure IT5|) that this differ- 
ence affects mostly the external radius. And we precisely 
note that that the radius of ui Cen is far too large, given 
its mass (Figure |2"2")) . 

Finally, the strong ellipticity of to Cen is likely the 
remnant of an initial difference from spherical symmetry, 
that has not enough time to vanish yet. 

9 CONCLUSIONS 

We are going to summarize the major results obtained in 
this work, and then list the directions in which it would 
be good to develop the research, in order to improve the 
theory and extend the range of its applications. 



9.1 Results obtained 

By making two hypotheses: (1) isotropy of the veloci- 
ties and (2) equal masses, we have obtained the system 
of equations (2.25), which, supplemented by the right 
boundary conditions, allows us to compute the evolution 
of a cluster from a given initial state. We have looked 
for an homologous solution, i.e. a model that remains 
self-similar, the evolution being limited to the scaling of 
the various physical quantities; we have found that this 
model exists and is unique. Its mass and its radius are 
finite. The projected density matches the observations as 
much as their precision allows. Near the center, the ho- 
mologous model yields unexpected properties: the central 
density and the central potential are infinite; a continu- 
ous flux of negative energy goes toward the center where 
it is absorbed by the formation of multiple stars. Fur- 
thermore, the escape of stars through the boundary of 
the cluster leads to a linear decrease of the total mass 
with time. 

Then, we have supposed that the main population of 
the cluster is mixed with a second population, much less 
numerous, of stars with a different mass, and we have 
computed the distribution of this secondary population 
and its escape rate. These results allow for an approxi- 
mate solution in the case of any mass function. By com- 
paring with the observations, we obtain in particular the 
approximate mass-luminosity relation of the stars in clus- 
ters. 

Finally we have studied the evolution of clusters close 
to the homologous model. We found that if the central 
density is finite, it rapidly increases and becomes infinite 
after a certain time. Generally, any difference with the 
homologous model decreases with time; a cluster becomes 
almost identical to the homologous model after the first 
third of its lifetime. In particular, these results lead to an 
estimate of the age of the globular clusters. 

9.2 Desirable improvements and extensions 

1. The most debatable hypotheses of the present model 
is indubitably the isotropy of the velocities. We 
could, as a first step, make the distribution function 
more general by adding an anisotropic term, small 
compared to the main term; this would allow us to 
quantify the error made when using the isotropy hy- 
pothesis. In a more ambitious second step, we would 
eliminate any restriction for the distribution func- 
tion f(E,A); but in this case, the equation of lo- 
cal evolution would become extremely complex (see 
Rosenbluth. MacDonald fc Juddll957l Equation 31). 



2. In the study of the approach to the homologous 
model made in Chapter VII, we assume that the 
cluster is already similar enough to the homologous 
model that we can linearize the equations of the dif- 
ferences. It would be very interesting to also study 
the evolution of clusters far from the homologous 
model; this would allow us to describe the beginning 



42 



M. Henon, translated by F. Renaud 



of the evolution of a cluster, from a given initial state. 
To do so, we would have to come back to the real 
[T.N;, physical] variables and directly solve (2.25), 
by computing step by step the successive forms taken 
by the cluster in timj^l. The results could be applied 
to a detailed explanation of the structure of ui Cen 
and to a precise estimate of its age. They would 
also allow for a more rigorous demonstration of the 
statement that clusters tend toward the homologous 
model, whatever their initial state is. 

3. The previous paragraph implies the knowledge of the 
initial state of the clusters, therefore the building of 
a more or less approximate theory for the formation 
of the clusters and the initial phase during which the 
steady state is established. We could assume, for 
example, that the protocluster is an homogeneous 
sphere with a density slightly higher than the critical 
density (therefore, favoring a contraction), with neg- 
ligible internal velocities. The evolution equations for 
this initial phase would be completely different from 
those we have considered here; on the one hand, the 
steady state is not reached and thus the distribution 
function does not only depend on E and A, and on 
the other hand, because this phase is very short, the 
effect of perturbations can be neglected. Fortunately 
this implies that there is no need for the knowledge 
of the time of formation of the stars. 

4. The hypothesis of equal mass stars is likely to be far 
from reality. In Chapter VI, we have obtained a more 
general solution but at the cost of a quite arbitrary 
approximation. Hence, we should extend the inves- 
tigations to the case of any mass function. But then 
a big difficulty appears: the escape rates of the stars 
of different masses not being the same, their relative 
importance always changes and the structure of the 
cluster varies too; as a consequence a homologous so- 
lution cannot exist. In other words, if one introduces 
the mass as an additional variable, one must also 
consider the time variable, and the complexity of the 
calculations is suddenly hugely increased. 

There are two possibilities to overcome this. First, we 
could seek a better model than the homologous one, 
that would allow us to separate the time variable; 
this model would necessary count, in addition to the 
two dimensional parameters, one or more other pa- 
rameters depending on time. Second, the mass func- 
tion could be initially described as a sum of [T.N;. 
Dirac's] S distributions. The simplest case of the sum 
of two 5 distributions (i.e. a cluster made of the mix 
of stars of two different masses, in any proportion), 
would already be a great improvement and would al- 
low us to estimate the actual effect of a spread of 
masses. 



5. We have only considered the value A = 1/3 of the 
power of the radius-mass relation (3.11), with the 
perspective of globular clusters. The calculations 
could be re-done with different values of A, and in 
particular, with A = oo, which correspond to the 
isolated cluster, and seems to better fit the cases of 
galaxies and galaxy clusters. The other values of A 
do not seem to correspond to existing objects, but 
could be useful to describe some intermediate states. 

6. The hypothesis of spherical symmetry is valid for 
most of the clusters; but some are clearly elliptical. 
Furthermore, we could foresee the application of the 
models to other objects (see point 8 below). Thus, 
we should study the non-spherical models. The el- 
lipticity can originate from a global rotation; a first 
approximation would consider that the rotation is 
small and only implies a correction term. We would 
likely find a model similar to the homologous one, 
but slightly flattened. It would be particularly inter- 
esting to find whether the escaping stars take away a 
little or a lot of angular momentum, and thus, if the 
flattening increases or decreases. 

We should also study the case of an angular momen- 
tum distribution without spherical symmetry, but 
still with no global rotation. Such an asymmetry 
would probably decrease, with a relaxation time of 
the same order as those found in Chapter VII for the 
radial differences. 

Finally, it would be good to consider the effect of the 
asymmetry in the galactic field (see Chapter III) . But 
this effect is likely not so important for the applica- 
tions, because it only affects the external regions, 
which are inv isible because of their very low density 
(|Kindll961ah . 



10 A similar calculation has been done, in the much simpler 
case o f a homogeneous plasma, by MacDonald, Roscnbl uth fc Wond 
|1957t) . 



7. From the observational side, it seems that much re- 
mains to be done. Accurate counts per luminosity 
class, as those made by Sandage for the cluster M3, 
would be very useful for other clusters. Furthermore, 
errors could probably be reduced by repeating the 
counts for se veral photographs of the same cluster. 
(lTavlerlll95l noted big differences between the plates 
obtained with two different telescopes.) 

The comparisons of accurate observation with a suf- 
ficiently elaborated theory would particularly allow 
for: 

• a precise physical determination of the mass- 
luminosity relation; to push this relation as far 
as possible toward the low luminosities, one 
should particularly study the closest clusters. 
In the future, space observatories would allow 
us to reach much higher magnitudes; 

• an estimate of the age of the clusters, and per- 
haps information about their formations, when 
analyzing the differences between the present- 
day structure of the clusters and their theoret- 
ical final structure. However, these differences 
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are small for most of the clusters, so that this 
analysis would require a high precision for the 
observations. One should start with the most 
massive clusters, where the largest differences 
remain; 

• a study of the galactic field and of the orbits of 
the globular clusters. 

Here, we have only focussed on globular clusters. The 
models obtained could be applied to other objects: 
galactic clusters [T.N.: star clusters in the disk of 
the galaxy, as opposed to the globulars that are fur- 
ther away in the halo], elliptical galaxies and cores 
[T.N.: central, spherical regions] of spiral galaxies, 
galaxy clusters. But in every case, difficulties appear, 
making the application of the theory more doubtful 
or more complex. That is, we think that the the- 
ory should first focus on describing properly the rel- 
atively simple case of globular clusters (which it is 
still far from doing) before extending its ambitions 
to more involved problems. 

The galactic clusters would be extremely interesting, 
for many reasons: being closer, less massive stars are 
visible there; their ages are not all the same and thus 
we can observe all the phases of the dynamical evo- 
lution; and it is possible to measure individual veloc- 
ities. Unfortunately, there are two major difficulties: 
the number of stars is small, so that statistical fluc- 
tuations forbid an accurate estimate of the projected 
de nsity and of t he other quantities; and above all, 



It seems worthless to underline the importance of the re- 
sults that could be obtained: the understanding of the 
collective dynamical phenomena is directly linked to the 
solving of the major cosmogonical and cosmological prob- 
lems, and thus to our understanding of the Universe. We 
hope that the present essay, although it has a very lim- 
ited scope, will contribute to the illumination of the new 
perspectives and give to others the wish to explore them. 
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as 



Spftzerl dl958h showed, the evolution of the clus- 



ter is largely influenced by passages near interstellar 
clouds. 



Finally, we note some minor issues: the enhancement 
of the computation of the perturbations between the 
stars; the study of the accumulation of central en- 
ergy as multiple stars and the final evolution of the 
cluster (see the end of Chapter V); the influence of 
the motion of the cluster within the galaxy on the es- 
cape of stars (this could be simplified as a three-body 
problem: galactic center, cluster, star). 



One sees that there is no shortage of work. Stellar 
dynamics is a rising science, where almost everything re- 
mains to be done. It looks today much more like a col- 
lection of isolated attempts than like a homogeneous doc- 
trine. These attempts have only led to isolated results, 
with limited range; the major problems, far from being 
solved, have remained almost untouched. 

As we have seen, the issues do not arise from a poor 
knowledge of the physical mechanisms involved, but only 
from the complexity of the calculations. But today, a new, 
extremely powerful, weapon exists to overcome this kind 
of issues: the electronic computers. That is, we think that 
stellar dynamics should now leave the relative state of ne- 
glect where it has been left and spark off the interest and 
the long-term efforts of a growing number of researchers. 



44 M. Henon, translated by F. Renaud 



List of the main notations 
The number is those of the equation where the notation is used first. 



a 


z.y 


77 

r 


2.26 


T 

L e 


0.0 


q 


O in 

z.iy 


t 


1 

z.l 


Zj 


K 1 
0.1 


A 

A 


2.2 


r 


z.ol 


L 


O OA 


Q 


z.z4 


1 


z.24 









Z.OO 




1 

0.1 


m 


1 

Z. 1 




z.oZ 


J. 


Z.6Z 


Q>i 


A Q 
4.0 


c 


9 ^fi 


r 2 


\J.O 


M 


9 41 


V 


9 7 


TJ 


^ .0 








6.6 


9 


2.12 


M 


2.43 


r e 


2.8 


u 


2.32 


13 


2.31 


c 


6.10 


H 


4.26 


M e 


3.12 




2.13 


U e 


3.5 


1 


2.31 


D 


2.24 


H 


4.31 


M P 


2.51 


R 


2.24 


U G 


3.1 


72 


6.5 


D 


2.32 


K 


4.33 


Mr 


2.45 


R 


2.32 


U 


7.1 


6 


6.22 


Dp 


2.47 


K 2 


7.10 


M 


2.38 


R e 


3.12 




5.7 


A 


3.11 


Dp 


2.48 


K D 


4.7 




3.3 


s 


7.23 


V 


2.3 




6.4 


E 


2.2 


L 


2.41 


n 


I 


S 


2.27 


X 


2.1 


P 


2.5 


E 


2.31 


L 


2.43 




6.24 


s 


2.32 


y 


2.1 


Pp 


2.46 


/ 


2.2 






P 


6.21 






z 


2.1 







The symbol ' indicates the derivation with respect to E or E. 
The subscript e corresponds to the values taken by the quantities at the boundary of the cluster. 
The symbols in capital letters (F, D, R, ...) generally represent the "normalized variables"; the bold font (F, D, R, ...) 

represent the "canonical variables" (see Chapter II) . 
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